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1.   INTRODUCTION 

The  principal  result  of  this  thesis  is  the  development  of  an 
algorithm  to  solve  the  large  sets  of  linear  equations  that  arise  in  the 
solution  of  elliptic  partial  differential  equations  by  implicit  numerical 
techniques.   The  algorithm  is  an  important  development  for  the  following 
reasons: 

1.  The  algorithm  is  always  convergent. 

2.  The  algorithm  requires  no  set  of  parameters. 

3.  The  algorithm  is  not  dependent  upon  special  properties  of 
the  partial  differential  equation  or  the  matrix  of 
coefficients  of  the  resulting  linear  equations. 

!+.   The  algorithm  has  proved  to  be  rapidly  convergent  on 
a  variety  of  test  problems.* 

5.   The  use  of  this  algorithm  in  conjunction  with  a 
factorization  due  to  Y.  J.  Kim  [11]  promises  to  be 
the  first  effective  iterative  method  to  solve  the 
systems  of  equations  generated  in  the  application  of 
finite  element  techniques. 

The  algorithm  developed  in  this  thesis  applied  in  conjunction  with  a 

factorization  is  the  only  existing  technique  with  all  of  these  properties. 


^Although  a  priori  estimates  of  the  rate  of  convergence  can  be  established, 
they  depend  upon  the  extreme  eigenvalues  of  the  iteration  matrix.   If 
accurate  approximations  of  these  eigenvalues  are  available  and  are  used  in 
the  algorithm,  then  the  algorithm  reduces  to  the  method  of  Dupont,  Kendall 
and  Rachford  [5]  with  optimal  parameters.   Thus  in  this  case  the  rate  of 
convergence  can  be  computed  as  described  in  [5].   The  application  of  the 
algorithm,  however,  assumes  no  knowledge  of  the  eigenvalues  of  the  iteration 
matrix.   Therefore  no  a  priori  bound  for  the  rate  of  convergence  is 
computable. 
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Since  the  application  of  direct  methods  to  the  large  systems  we 
are  considering  is  either  impossible  or  impractical,  iterative  methods  are 
usually  used.   Unfortunately  the  standard  iterative  techniques  are  also 
deficient.   The  Jacobi  method  is  quite  slow  even  on  some  simple  problems. 
The  Gauss-Seidel  [19]  method  can  be  shown  to  be  simply  a  factor  of  two 
better  than  the  Jacobi  iteration  regardless  of  the  number  of  equations. 
Successive-Over-Relaxation  [19 ]  is  dependent  on  the  problem  and  may  become 
quite  slow. 

The  Alternating  Direction  Implicit  iteration  [20]  is  more  efficient 
than  either  Gauss-Seidel  or  Successive-Over-Relaxation.   However,  it  requires 
the  selection  of  a  set  of  iteration  parameters  to  be  applied  cyclically 
during  the  iteration.   In  general,  there  is  no  theoretical  basis  for 
selecting  a  set  that  will  give  rapid  convergence.   Moreover,  the  method 
loses  its  fast  convergence  when  the  region  of  definition  of  the  boundary 
value  problem  is  irregular. 

The  algorithm  developed  in  this  thesis  is  based  on  factorization 
techniques.   To  elucidate  the  difference  between  this  algorithm  and  other 
algorithms  based  on  factorizations,  we  shall  give  a  brief  explanation  of 
factorization  methods. 

The  solution  of  a  system  of  equations  Ax  =  q  by  factorization 
techniques  can  be  separated  into  two  segments.   The  first  step  is  the 
definition  of  an  auxiliary  matrix,  B,  which  is  chosen  so  that  A+B  =  LU 
where  L  and  U  are  sparse  lower  and  upper  triangular  matrices  respectively. 
In  practical  computer  programs  the  elements  of  B  are  not  computed;  instead 
the  elements  of  L  and  U  are  computed.   Then  by  using  the  factors  L  and  U, 
systems  of  equations  whose  matrix  of  coefficients  is  A+B  can  be  easily  solved. 
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The  second  part  of  the  method  consists  of  an  algorithm  which  uses 
an  iteration  of  the  form 

(A+B)xn+1  =  (A+B)xn-  (Axn-q),  n  =   0,  1,  2,    ... 

to  generate  a  sequence  of  vectors  {x  }  which  converges  to  the  solution  of 

Ax  =  q.   The  convergence  of  the  iteration  is  dependent  upon  the  definition 

of  the  auxiliary  matrix  (through  the  definition  of  L  and  U)  and  a  sequence 

of  parameters. 

One  of  two  types  of  parameters  is  employed  in  the  iteration.   The 

first  type  may  be  used  to  alter  the  auxiliary  matrix  for  each  n.   There  has 

been  no  analysis  of  the  resulting  iteration, 

(A+B  )x  ,  =  (A+B  )x   -  (Ax  -q). 
n  n+1       n  n      n 

Therefore  there  is  no  theoretical  basis  for  the  selection  of  such  a  set  of 

parameters. 

The  second  set  of  parameters,  denoted  by  x  ,  multiplies  the  term 

(Ax  -q)  in  the  iteration  as  follows 
n 

(A+B)x  _  =  (A+B)x  -x  (Ax  -q). 
'  n+1        n  n   n  ^ 

In  previous  algorithms  the  definition  of  the  parameters  x  is  based  upon  a 

%Ay  x^                                *\Ax  x  ^ 
lower  bound  for  min  yjrrhn ^  and  an  upper  bound  for  Max  prrfrrl <>  I 

x^0  <{A+B)x,x>  **  x£0  <{A+B)x^x>  ' 

equivalently  they  are  based  on  lower  and  upper  bounds  for  the  extreme 
eigenvalues  of  (A+B)  A.   These  bounds  must  be  calculated  every  time  the 
matrix  of  coefficients  A,  is  altered.   Furthermore,  the  rates  of  convergence 
of  the  previous  algorithms  are  dependent  upon  the  accuracy  of  the  bounds. 
The  new  algorithm  requires  no   knowledge  of  the  eigenvalues  of  (A+B)  A. 
Moreover  the  rate  of  convergence  is  improved  by  quantities  calculated  during 
the  execution  of  the  algorithm  in  such  a  way  that  the  rate  of  convergence 
approaches  the  optimal  rate. 
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The  new  algorithm  is  also  superior  to  algorithms  which  vary 
the  auxiliary  matrix  for  each  iteration.   Its  primary  advantage  is  that 
there  is  a  mathematical  basis  which  guarantees  convergence.   A  second 
advantage  is  that  only  one  auxiliary  matrix  is  used.   Thus  only  one  pair 
of  factors  L  and  U  is  defined  and  consequently  less  work  is  required. 

Although  no  knowledge  of  the  eigenvalues  of  (A+B)~  A  is  necessary, 
information  about  the  eigenvalues  is  useful.   Furthermore,  the  rate  of 
convergence  is  governed  by  the  extreme  eigenvalues.   Therefore  the  selection 
of  the  auxiliary  matrix  is  still  important.   For  matrices  which  arise 
in  the  solution  of  elliptic  partial  differential  equations  by  finite 
difference  methods,  the  symmetric  auxiliary  matrix  defined  by  Stone  [15] 
(see  section  2.2)  seems  to  be  the  best  choice  of  the  known  factorizations. 

We  shall  now  summarize  the  results  of  this  thesis. 

In  Chapter  3  it  is  shown  that  for  arbitrary  A  and  A+B  the  iteration 

(A+B)x  , .  =  (A+B)x  -t  (Ax  -q) 
n+1        n     n       n  ^ 

is  convergent  with  t  a  constant  t  iff  (if  and  only  if)  the  real  part  of  each  of 

-1  t 

the  eigenvalues  of  A  B  is  greater  than  p  -  1.   By  modifying  this  result  for 

positive  definite  A  and  Hermitian  B,  it  is  shown  that  the  iteration  is 

convergent  with  t  =1  iff  A+2B  is  positive  definite.   This  strengthens  the 

long  established  result  that  when  A  is  positive  definite  and  B  is  Hermitian, 

the  iteration  is  convergent  with  %     =   1  if  A+2B  is  positive  definite,  see 

Wachspress  [20].   Next  it  is  shown  that  if  A  and  A+B  are  positive  definite 

then  the  iteration  is  convergent  with  a  constant  value  t  iff 

2 


0  <  T  < 


Max  X. ((A+B)_1A) 


where  \  ((A+B)"  A)  is  an  eigenvalue  of  (A+B)-1A. 
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In  section  3.3  the  analysis  by  Dupont,  Kendall  and  Rachford  [5] 

concerning  the  selection  of  a  best  sequence  of  x  ' s  to  use  when  the 

eigenvalues  of  (A+B)  A  are  real  is  presented.   The  results  are  extended 

in  section  3-h   establishing  a  best  sequence  of  x  's  to  use  when  the 

n 

eigenvalues  of  (A+B)  A  are  complex. 

In  Chapter  k   an  algorithm  based  on  the  iteration 

(A+B)x   .  *  (A+B)x  -x  (Ax  -q) 
n+1        n  n   n  *' 

is  developed. 

In  Chapter  5  an  approximation  of  the  smallest  eigenvalue  of 

(A+BQ)~  A,  where  B  is  the  symmetric  auxiliary  matrix  defined  by  Stone ,    is 

b  b 

derived  for  some  special  cases.   This  approximation  can  be  used  in  starting 
the  algorithm  of  Chapter  >+  when  BQ  is  the  auxiliary  matrix. 

b 


2.   PEELIMINARIES. 

Notation:  Throughout  this  thesis,  X. (M)  will  denote  an  eigenvalue  of 

the  matrix  M,  X   .    (M)  will  denote  the  eigenvalue  of  minimum  algebraic  value, 

\^     (M)  will  denote  the  eigenvalue  of  largest  algebraic  value,  and  p(M) 

will  be  used  to  denote  the  spectral  radius  of  M,  i.e.,  p(M)  =  Max  \.  (M)  . 

i   -1   I 

All  matrices  used  in  this  thesis,  even  those  specified  as  general, 
are  real  matrices. 

2.1  Origin  Of  The  Problem.   Consider  the  Dirchlet  problem 

j 

-  ^  (ai(x)  1^  -  ox;  (a2(x)  si:)  =  q(x)>  x  eD'  u(x)  =  {x)>  for(2.f) D 

where  x  =  (x  ,xp),D  is  the  interior  of  a  compact  region  with  boundary  c^D 
and  where  the  differential  operator  in  (2.1)  is  elliptic,  that  is  a, (x) 
and  ap(x)  are  strictly  positive  on  D.  A  common  method  of  solving  (2.1)  is 
to  approximate  the  differential  operator  with  a  difference  operator.  A 
system  of  linear  equations  results  and  the  problem  is  then  to  find  a  solution 
to  the  system  of  equations. 

The  system  is  derived  as  follows.   Assume  D  lies  in  the  first 
quadrant;  cover  D  with  two  families  of  parallel  lines, 

x1  =  jh      j  =  1,  .  .  . ,  n± 
and 

x2  =  kh      k  =  1,  . . .  ,  np  . 

Denote  the  points  of  intersection  (jh,  kh)  by  (j,k).   The  points  (j,k)  are 
called  the  mesh  points  or  gridpoints.   Let  bj),    be  the  polygon  formed  by 
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joining  those  mesh  points  (j,k)  for  which  one  of  (j+l,k)  or  (j,k±l)  is 
not  in  D.   Similarly  let  D  be  the  open  and  connected  set  hounded  by  dD  . 
At  each  gridpoint  of  D  the  differential  operator  is  approximated  by  a 
difference  operator  and  the  equations  are  in  one-to-one  correspondence  with 
the  gridpoints. 

For  simplicity  of  notation  we  shall  consider  the  region  D  to  be 
the  unit  square  bounded  by  the  x-.  and  xp  axes  and  the  lines  x  =1  and  xp=l. 
Let  x  =  (x1}1,   r^g,  ...,   x^,  x2 A,    ...,   xn^)T  where  n±   =   n?_  =  n.   We 
obtain  the  matrix  equation 

Ax  =  q  (2.2) 

defined  by  (Ax) .    *  q.  .  ,   where  the  subscript  corresponds  to  gridpoint 
(j^k)  and  q    is  derived  from  q(jh,kh)  and  the  boundary  conditions. 

The  difference  operator  can  be  chosen  so  that  matrix  A  is  diagonally 
dominant  and  positive  definite  with  a...  >  0  and  a. ,  .  <  0  for  i  /  j.   One 

method  for  deriving  such  an  A  is  to  approximate  3 —  (a.  (x)  3 —  u(x))  by 

i  i 

h"  {ai(x+  g-e  )   [u(x)-u(x+hei)  ]  +  a±(x  -  ^e±)      [u(x)-u(x-hei)  ]  } 

where  e  and  e  are  unit  vectors  along  the  xn  and  xp  axes  respectively.   This 
defines  A  by 

(Ax).,  n  B.  ,  X.  ,  ,  +  D.  .  x.  ,  ,  +  E  .  .  +  D. ,,  .  x.,_  ,  +  B.  ,  .x.  ,,.  , 

yj,k    j,k  j,k-l    j,k  j-l,k    j,k    o+l, k  j+l,k    j,k+l  j,k+l 

j^k  ss  1,  , . .  t   n,   where 
B •  k  =  -a2(jh,(k-  g-)h), 

Dj  k  =  -ax((j-  §)h,kh),  (2.5) 


and 

E.   =  a2(jh,(k-  |)h)  +  a2(jh,(k  +  |)h  +  a^j-  §)h,kh) 

+  ax{U  +  |)h,kh). 
When  1  =  1  or  k=l  we  define  D.  .  =  0  or  B.   =  0  respectively. 

The  following  properties  of  A  and  the  quantities  defined  in  (2.3)  will  be 
used  below. 


^,k  -Ed,K*a,k  +  Bj,kxj,*.i  +  Vxd-i,*  +  ViaV,*  +  Bj,k+ixj,k+i  , 


and 

E.   >  -2(B.  .  +  D   )  with  equality  holding  if  B.  ,  ^  0  and  D.  .  /  0 
J^k  ~~     J^k    j,k  J^-K         J^-K 

for  j=l,  2,    . . .  ,   n  and  k=l,  2,    . . .  ,  n. 

The  subscript  («j,k)  refers  to  the  gridpoint  (jh,kh). 

In  this  thesis  we  shall  study  the  solution  of  equations 

Ax  =  q 

by  factorization  techniques.  As  sescribed  in  the  introduction,  the  idea 

of  a  factorization  method  is  to  add  an  auxiliary  matrix,  B,  to  A  so  that 

A+B  factors  into  known  sparse  matrices  L  and  U.   An  algorithm  based  on  the 

iteration 

(A+B)x  .  =  (A+B)x  -t  (Ax  -q)  (2.5) 

n+1      '  11  n   n 

is  then  used  to  compute  a  sequence  of  vectors  {x  )  which  converges  to  the 
solution  x.   The  parameters  1     are  defined  by  the  algorithm. 
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A  factorization  is  defined  by  the  auxiliary  matrix  or  the  factors 
L  and  U.   We  shall  give  particular  attention  to  the  symmetric  factorization 
due  to  Stone  [15  ].   This  factorization  which  can  be  applied  to  any  matrix  A 
with  the  properties  listed  in  {2.1+) ,   is  defined  as  follows. 


where 


j,k  '   j,kXj,k-l   cj,kxj-l,k  H   j,kxj,k  ' 

(Ux)  .  ,  a  x .  ,  +  e.  ,x.r,  ,  +  t.   ,X.  ,  -,  > 

yj,k    j,k    j,k  j+l,k    j,k  j,k+l 


b .  ,  =  B .  ,  -ac .  .  _  f .  .  ,  ,  , 
j,k    j,k   j,k-l  j-l,k-l 


(2.6) 


c  .  1  =  D .  ,  -ab     e 
J,k    j,k   j-l,k  j-l,k-l 


j,k  j,k  d,k-l  j,k  j-l,k  '   j,k   j,k-l  J-l,k-l 

+ab     e 

j-l,k  j-l,k-l  , 


d.  e  .   s  D.,n  .  -ab  .  .  e  .  .  .  , 
j,k  j,k    j+l,k   j,k  j,k-l 


(2.7a) 


and 

j,k  j,k  ""  j,k+l"  °j,k  j-l,k 

The  quantities  B.  .  ,  D.  ,  and  E.  .  are  the  elements  of  A  defined  in  (2.4), 
j,k    j,k      j,k  y' 

and  the  parameter  a  is  used  to  vary  the  auxiliary  matrix  if  Stone's  algorithm 
[16]  based  on  the  iteration  (2.5)  is  used.   We  are  interested  in  algorithms 
that  don 't  vary  the  auxiliary  matrix  and  will  assume  a  =  1. 

In  addition  to  the  quantities  defined  in  (2.7a)  we  define 

Ci  k  "  bi  kei  k  1     and     Gi  k  =   Ci  kfi  1  k  '  (2'Tb) 

Since  there  should  be  no  confusion  as  to  whether  the  auxiliary  matrix 
is  defined  by  (2.6)  and  (2.7)  or  is  arbitrary,,  B  will  be  used  to  denote  either. 
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With  L  and  U  of  the  form  (2.6),  A+B  is  given  by 

((a+bW^  =  ((MWj,*  -  V^A-i  +  tJ,kej,k.ixj+i,k-i 

+  cj,kxj-l,k  +   j,k  H   3,k  j,k-l  "*  cj,kej-l,k',xj,k 

+  dj,k  ej,kxj+l,k  +  cj,kBj-l,kxj-l,k+l  +  dj,kfj,kxj,k+l  * 

(2.8) 

2.2  Known  Results   The  following  definitions  are  standard. 


T  \T 

Definition  1.1:  If  x=(x1,  . . .  ,  X  )  and  y~(y,,  ...,  y  )  are  n-component 

n     *       * 

vectors  then  ^c,y^  =  T\    x.y.  where  y.  denotes  the  complex  conjugate  of 

L=l  X  X       X 
y.  .  -NX,y>  is  called  the  innerproduct  of  x  and  y. 

jL, 

2 

Definition  1.2:      ||x||,   the  L2  norm  of  x,   is  defined  by  ||x||    =  ||x  ||  _      =  <x,x>     . 

2 

If  A  is  positive  definite  then  the  A- norm  of  x,  ||x||.,  is  defined  as<Ax,x>2. 


Dupont,  Kendall  and  Rachford  [Jj,]  have  proved  the  following  lemma. 


Lemma  2.1:   Let  A  and  A+B  be  positive  definite.   Suppose  there  exist  positive 
numbers  e  and  ep  such  that  for  all  nonzero  x  we  have 

<(A+B)x,x>  e  [el'e2]  (2'9) 

Then  for  0  <t  <2e  ~  the  sequence  {x  }  defined  by  (2.5)  converges  to  the 
solution  of  Ax=q.   Furthermore,  for  i  =  2{e     +   e  )~  the  A-norm  of  the  error, 
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e   =  x-x  ,  is  reduced  by  a  factor  of  at  least  (e  -e  )/(e  +e  )  in  each 
iteration. 

It  is  also  suggested  [k]   that  a  sequence  of  t   derived  from  a 
classical  Tchebychev  analysis  be  used.   Iteration  (2.5)  is  then  performed 
in  practice  by  using  a  technique  due  to  Stiefel  [k].      This  technique  will  be 
examined  in  Chapter  3  an(i  a  modification  of  it  is  used  in  Chapter  1+  in 
developing  a  new  algorithm. 

Bracha  [1]  has  derived  the  following  important  properties  of  the 
elements  of  the  factors  L  and  U  defined  in  (2.6)  and  (2.7). 

M   V<B.A<0, 

M      e^<D.>k<0, 

M      d.i,*>0> 

(a)  -1<  e   <  0, 

(e)  -1<  f    <  0, 

J  >& 

and 

(f)  1+e.  ,+f..  .  >  0. 

Bracha  has  also  shown  in  [1]  that  for  the  factorization  defined  in 
(2.6)  and  (2.7)  A+B  is  positive  definite  when  A  satisfies  (2.4).   A  more 
concise  proof  is  included  here. 

Lemma  2.2:   If  A  satisfies  (2.k)   and  B  is  defined  by  (2.6)  and  (2.7)  then 
A+B  is  positive  definite. 


Proof:  From  (2.7)  it  is  seen  that  b  .  ,  =d..  ,  f .  ,  n  and  c.n  =  d.  -  .  e.  -  .  . 
v  "  j,k    j,k-l  j,k-l     j,k    j-l,k  j-l,k 

Thus  (2.6)  implies  (L)    =  d  •  •,  •  (u)     where  row  r  of  A+B  corresponds  to 

r,s    j  ,k.  s,r 
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T 
the  gridpoint  (j,k)  and  s  =1,  . . .  ,   r.  We  have  L  =  U  diag(d  ),   where 

j  )£■ 

T 

diagfd.  .  )  is  the  diagonal  matrix  with  d.  ,    along  the  main  diagonal  and  U 
J,k  ot* 

is  the  transpose  of  U.      Using  property  (c)   from  above,  we  have 

A+B  =  LU  =  UTdiag(d1/2)   diag(d1/2)U 
which  is  clearly  positive  definite. 


In  addition  Bracha  proved  the  following  theorem   [2]. 


Theorem  2.1;     Let  A  and  A+B  he  the  matrices  defined  in  (1.6)   and  (l.  10) 
respectively.      Let  e.       and  f         be  defined  hy   (1.10).      Let 


p,j,k =  -e. ;  -f. ;     J**  - x>  z>  - •• >  n> 

J,k        j,k 


and 


p  =  min  p         . 

Then 

<Ax,x>  ^  1 

Max      ttttsR — r     <      r* 

x?fo       <(A+B)x,a>       -         1  _  J^ 


««j  •  "^Ax.x)"  .  1 

and  min       yl .   *S       g        >     

x/0      <(A+B)x,x>        "        l  +  ki 


where  k  is  a  positive  constant. 


3-  CONVERGENCE  CRITERIA 


In  this  chapter  necessary  and  sufficient  conditions  are  derived 

for  the  convergence  of  (A+B)x  ,,=  (A+B)x  -t  (Ax  -q).   The  conditions  we 

n+1       n  n   n 

obtain  do  not  seem  to  he  know  in  the  literature.   The  conditions  are  first 

developed  for  arbitrary  A  and  B.   For  this  case  it  is  shown  that  the 

iteration  is  convergent  with  a  constant  x  iff  i   is  selected  from  a  certain 

region  in  the  complex  plane.   This  result  is  modified  for  the  case  when  A 

is  positive  definite  and  further  modifed  for  the  case  when  both  A  and  A+B 

are  positive  definite.   When  A  and  A+B  are  positive  definite,  a  best  constant 

value  of  r  is  established. 

Section  ^.k   establishes  a  best  sequence  of  t  's  to  use  when  the 

eigenvalues  of  (A+B)"  A  are  complex.   The  derivation  follows  that  of  the 

known  results  presented  in  Section  3-3;>  which  establish  a  best  sequence 

of  t  's  when  the  eigenvalues  of  (A+B)  A  are  real. 
n 

3.1  Convergence  Criterion  When  A  and  A+B  Are  Arbitrary 


In  the  following  two  lemmas  convergence  criteria  are  derived  in 

terms  of  t  for  (A+B)x   =(A+B)x  -t(Ax  -q).  (3-1) 

n+1      n    n 


Lemma  3.1:     Iteration  (3.1)  is  convergent  iff  Re(X. (A  B))  >  _  -  1. 


th 
Proof:     Let  x  be  the  exact  solution  of  Ax=q.   Let  the  error  in  the  n 

iterate  be  given  by  e   =  x-x  .   Subtracting  (3-l)  from 

n      n 

(A+B)x=  (A+B)x-T(Ax-q) 

yields    (A+B)  (x-Xn+1)  =  (A+B)  (x-XQ)-T(A(x-xn)) 

or  e   =  (A+B)"  (A+B-rA)e   . 

n+1  n 
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Ik 

The  matrix  (A+B)"  (A+B-tA)  is  the  iteration  matrix  of  the  iteration  (3.1)- 
It  is  well  known  that  the  iteration  is  convergent  iff 

p(  (A+B)"1(A+B-tA))<1.  (3.2) 

The  iteration  matrix  can  he  rewritten  as 

(A+B)"1  (A+B- xA)  =  (I-t(A+B)"1A) 

=  (I-t(A-1(A+B))-1)  =  (I-Td+A^B)-1), 
where  I  is  the  identity  matrix. 
Thus  (3.2)  can  he  rewritten  as 

^.(I-Td+A^B)"1)  eS  (0,0;  1) 


where 


5(u>v;R)  =  {(u,v):(u-x)  +  (v-y)  <R  } 


Hence  (3.1)  is  convergent  iff 

^(-Ll-d+A^B)"1)  e  5  (0,0;  4)  >  (5-3) 

^((I+A^B)"1)  e  5  (-i,0;  ±)    , 

or  \. (l+A~  B)  has  real  part  greater  than  x/2.   It  now  follows  immediately 
that  (3«l)  is  convergent  iff 

Re(\  (A_1B))>  1-    -1.  (3. k) 

Corollary  3. 1.1:  If  A  is  positive  definite  and  B  is  Hermitian,  then  (3-1) 
is  convergent  iff  A+B-  ^-A  is  positive  definite. 


Proof:   If  A  and  B  are  Hermitian,  each  \.(A~  B)  is  real  and  (3.4)  becomes 

?,.(A_1B)>  -L-l  . 
This  is  equivalent  to  K.  (l(l-  i)+A  B  )X>  ■   But 

1(1-  1.  )+A_1B  =  A"1(A(1-  1  )  +B)  =  A_1(A+B-  -£-A). 
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Thus  (3. if)  is  equivalent  to  X.  (A~  (A+B-  -*-A))  >  0.  (3.5) 

Since  A-1(A+B-  ^-A)  is  similar  to  A"1/2  (A+B ^-A)A_1/2   , 

(3.5)  may  "be  restated  as 

K  A  '   (A+B-  -*-A)A~    x,x)>>0,  for  all  nonzero  vectors  x. 
Therefore 

<  (A+B J-A)A"  '  x,A~  '   x  >  >  0,  for  all  nonzero  x. 

Letting  y=A~    x  we  obtain  K   (A+B-  -—-A)y,y  ^  >  0,  for  all  nonzero  y. 

completes  the  proof. 


Hence  (3.  l)  is  convergent  iff  A+B-  -^A  is  positive  definite,  and  this 


For  x  =  1  the  corollary  implies  (3«l)  converges  iff  A+2B  is 
positive  definite.   This  strenghtens  the  previously  known  fact  that  a 
sufficient  condition  for  the  convergence  of  (3.1)  is  that  A+2B  he  positive 
definite. 

3.2  The  Positive  Definite  Case 

We  now  modify  the  above  results  for  an  A+B  which  is  positive 
definite,  as  it  is  in  the  factorization  due  to  Dupont,  Kendall  and 
Rachford  [ 5]  and  the  symmetric  factorization  due  to  Stone  [ 15] . 


Lemma  3.2:     If  A  and  A+B  are  positive  definite  then  (3-1)  is  convergent 

2 
iff  0  <  x  <  . 

Proof:         (3«l)  is  convergent  iff 

I  ?,.(I-t(A+B)_1A)  I  <  1. 
Since  \.((A+B)~  A) >  0  ,   this  condition  is  equivalent  to 

1"tNiax((a+b)"1a)>  -1 

2 

\  ,    and  this  completes  the  proof. 

3--J-A    1 


W(A+B  A> 
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We  can  make  a  further  refinement  by  deriving  a  "best  value  for  x 
as  seen  in  the  next  lemma. 


Lemma  3.3:     If  A  and  A+B  are  positive  definite  then 
2 


X,   = 


*  "    Am((A+B)-1A)+Xmln((A+B)-1A) 
is  the  best  choice  for  t  in  (3.1). 


-1, 


Proof:     The  iteration  matrix  of  (3.1)  is  T  =  (I-t(A+B)~  A).   Thus,  since 
all  the  eigenvalues  of  (A+B)~  A  are  positive,  the  eigenvalues  of  T  are 
decreasing  as  a  function  of  x.   Clearly  the  spectral  radius  of  T  is  minimize' 

with  respect  to  1,   when  |  ^-ji-C^)  I  =  I  ^mAX^t^  I  '   ±'e' > 

I  l-T\./rA..((A+B)'1A)  I  =  1-t\  .  ((A+B)_1A)  , 
1     MAX.     '     '       mm  ' 

or       2  =  t[  \MAX((A+B)-1A)+\min((A+B)"1A)]. 

Therefore  t,  =  2/[\,,/.„((A+B)"  A)  ]  is  the  best  choice  of  x  as  claimed. 
b      m  2X    .  ((A+B)"^) 


Moreover  for  x  =  x  ,p(T  )  =  1- 


min 


1  ,'",'"r"      T^.MA+Byh^X   .  ((A+B)"^) 

MAX  v   '        '     mm 


There  is  a  strong  relation  between  Lemmas  3. 2  and  3-3  an(i  Lemma 
2.1.   Lemma  3.2  states  a  necessary  and  sufficient  condition  while  Lemma 
2.1  states  only  a  sufficient  condition.   The  hypothesis  of  Lemma  2.1, 
however,  includes  a  supposition  that  positive  numbers  e  and  e  exist  such 
that  for  all  nonzero  x 

<  Ax,x  >        e  [e,  ,e  ]  . 
<  (A+B)x,x  >  X     d 

This  is  certanily    valid  when  both      A  and  A+B  are  positive  definite. 

In  fact 


e,    = 


K»-   '(A) 


mm 


1         W  (A+B) 
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and 

XMAX  (A) 

e2    ~     [a+bJ 

iran 
are  sufficient.   The  next  lemma  shows  precisely  the  connection  of  Lemmas 

2.1  and  J. 2. 

Lemma  3. k:        Let     A  and  A+B  be  positive   definite.      Then 

<    AX.X     >  r  n 

<  (A+B)x,x  >  €    Lel;62J  for  all  nonzero  x 

iff  \.((A+B_1A))   e    [e±,e2]   . 
Proof:  The   conclusion  may  be  restated  as 

1  <  ay,y  > 

X   .    ((A+B)     A)    =    Min         j   ,.       . <—  .     c. 

minvv  ,  <   (A+B)y,y  >  (3.6) 

and  y  v 

2_  <  ay,y  > 

\        ((A+B)'  A)    -  Max       . 

™  y^o  <  (A+B)y,y  > 

Now  since      (A+B)~  A    is   similar  to  A  •'<    (A+B)     A  '^     which  is  Hermitian 

(3.7) 

-1a ^    _,         ^V2    rAo.^-l4V2    >    _     m,„  <A1/2    (A+B)-1Al/2> 


\_.  m  ( (A+B)     A)    -  \_.  „  (A  '       (A+B)     A  '       )    =     Min 

x/o 


mm  man 


x 


.  1/2    ,        ,_i  1/2    -  -1/2  1/2 
Letting  x   =    LA  (A+B)     A  J  Ay  and  using  "the   selfadjointness 

of    [A1/2    (A+B)"^2    j"1/2  and  A1/2      (3-7)  becomes  \   .    ((A+B)"^)      = 

min 

Mnn      <  Ay,y  > 
$0   <  (A+B)y,y  >  ' 

Similarly  ^((A+B)'^)  =  Max  <   (A+B^y  >      This  completes 
the  proof. 

3.3  The  Use  Of  A  Sequence  Of  t's  When  A+B  Is  Positive  Definite 

There  are  no  new  results  presented  in  this  section.   However, 
Section  J.k   and  Chapter  k   are  developed  from  the  results  described  here. 
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In  this  section  a  "best"  sequence  of  t  's  to  be  used  in  the 
iteration  (A+B)x  _  =  (A+B)x  -Tn(Axn-q)  ,  is  derived  for  the  case 
when  both  A  and  A+B  are  positive  definite.   A  "best"  sequence  is  one  that 
minimizes  the  spectral  radius  of  the  iteration  matrix. 

As  in  Lemma  3.1  the  error  vector  e  is  multiplied  at  each  step 

by  (I-t  (A+B)"  A).  Hence 

N-l  , 

6N  =nio   ^-Tn(A+B)'  A)  e0  '  (5'8) 

To  make  the  iteration  converge  as  rapidly  as  possible,  the  t  's  should  be 

chosen  to  minimize  the  spectral  radius  of 

N-l  1 

"n  ■  nSo  <WA+B)"  A)- 

The  eigenvalues  of  M^.  are 

N-l  1 

W  -&>  [I-Vi(<A+B»"  A'! 

and  since     \.((A+B)~  A)    is  real,   •  \.(^L)   is  real.   The  problem  of 
minimizing  the  \.  (HL  )  thus  reduces  to  selecting  a  sequence  of  t  's  that 

minimizes  the  maximum  absolute  value  of 

N-l 
PN(X)  =n&  (l"TnX)  (5'9) 

on  the  interval  X   .  ((A+B)"  A)  ^  X  ^  \,,._-((A+B)"  A).   No  normalization 

mm     '        J  MAX 

is  necessary  since  it  is  clear  that  all  P„  of  the  form  (3.9)  satisfy  P  (o)=l. 

By  classical  Tchebychev  analysis  it  is  seen  that  the  desired 
polynomial  is  given  by 

p  (y)      _  V  (?W+:SlAX  -2X)/(XMAX-Xmin))   ,  (3.IO) 

Nv  mm  MAX  '     TdAX  mm' 


where  T  is  the  Tchebychev  polynomial  of  degree  N  and  \   .   and  K^jr  ^en0^e 
\da((A«)"1A>ana  ^((A+l 
of  PH(X)  on  ft^,^]  is 


^•min((A4B)"  A)  and  ?\^^((A+B)~  A)  respectively.   The  maximum  absolute  value 
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|~         /   Nnin  +XMAX  \     ~~| 

L         \     MAX       mm/     J     .  (3-H) 

Note  that  if  \  .  ^  0  <  \.flY  then  the  maximum  of  P  must  he 
greater  than  or  equal  to  1,  since  P  (0)=1.   However  a  and  A+B  positive 
definite  implies  that  X^  and  \^^  are  positive  and  therefore  assures 
the  reduction  of  the  error  vector  when  {t  }  is  defined  by  (3.9)  and  (3.10) 

There  are  two  ways  of  implementing  a  sequence  of  x  '  s  so  that 
(3.9)  is  the  polynomial  given  in  (3. 10).   In  the  first  method  an  N  is 
selected  and  the  sequence  x  n=0,l. . .  ,N-1  is  computed  as  the  reciprocal 

of  the  roots  of  T   /  min  MAX   \  .   The  sequence  is  then  used  repeatedly 

^MAX'Nnin 

until  the  norm  of  the  residual  vector,  r  =  Ax  -q,  is  sufficiently  small. 

'   n     n  J 

A  sequence  of  x's  derived  as  ahove  is  called  a  Tchebychev  sequence  of 
length  N. 

The  other  method  of  implementing  the  sequence,  which  was  first 
suggested  by  Stiefel[l|],  uses  the  recursive  property  of  the  Tchebychev 
polynomials.  If  this  technique  is  used  then  every  x  is  the  result  of  a 
Tchebychev  sequence  of  length  n.  if  the  first  method  is  used  then  x  is 
the  result  of  a  Tchebychev  sequence  of  length  n  only  when  n=N. 

The  Stiefel  sequence  is  defined  as  follows.   Suppose  x  is  the 
result  of  using  a  Tchebychev  sequence  of  length  k.   Define  x    as 

where 

(\lAX  "  WW*'  \^>  "  ' 


with 


2       (A+B)"1  (q-AxJ 
5x  =- ~ V    ;    V4   0y 

0   Nnin  +  ""MAX 
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for 


Y  _  HlAX  +Xmin 


*"MAX  ~^min 
then  x.    is  the  result  of  using  a  Tchebychev  sequence  of  length  k+1. 

K+l 

In  general  \  .   and  \.,.„  are  not  known.  Tchebychev  minimization 
°       mm     MAX 

is  typically  implemented  by  finding  an  interval  (e  ,eg)  containing 
(>,min,>mx)  and  using  e±   and  e2  instead  of  X^  and  k^  respectively. 

In  place  of  (3»9)  we  have, 


e.+e^-2X 


N-l 

p„(x)  =  n  _  (i-T  x)  = 


And  the  Stiefel  sequence  is  defined  by 


H 

e2"el 

J 

t  r 

6l+e2 

1 

JM 

e„-e„ 

Vi  =  \  +  5 \  > 

Wk(y)  Tk  1(y) 

5  *k  -  (e2-ei)Tk+1(y)   (A+B)   (^V  +  T^y)   Vl  , 


S  x^  = 


(A+B)_1(q-Axn) 


0    e!+e2  °* 

where  y  =  (e  +e  )/(e2-e  )   . 

3.1+  Obtaining   A  Sequence  Of  t's  When  A+B  Is  Not  Positive  Definite 

In  this  section  we  shall  derive  a  best  sequence  of  t's  to  be 


used  in 


(A+B)x  _  =   (A+B)x  -t  (Ax  -q) 
n+1      '   n  n   n  ^ 
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-1, 


when  (A+B)  A  has  complex  eigenvalues.  Before  proceeding,  it  is  necessary 
to  define  what  "best"  means.   The  iteration 


(A+B)x  , ,  -  (A+B)x  -t  (Ax  -q) 
n+1        n  n   n 


N-l 

n=  0 


is  most  rapidly  convergent  when  p(r|IIn(l-T  (A+B)~  A))   is  minimized. 


This  is  accomplished  when 


N-l 


-1, 


Max  n   (1-t  X. (A+B)  A) 
'   n  i 
i  n=  0 


(3-13) 


is  minimized.   Since  not  all  of  the  X.    are  known,  the  t  could  he  selected 

1  n 

to  minimize 


Max 
zeR 


Vz> 


N-l 
where  P,T(z)  =  II   (l-x  z) 
N    a».0 


b-ik) 


and  R  is  some  region  containing  all  of  the  eigenvalues  of  (A+B 
This  concept  is  formalized  with  the  following  definition. 


-1 


A. 


Definition  3- !•  A  sequence  of  parameters  1     n  =  0,1,«««,N-1  used  in 
iteration  (2.5)  is  called  a  best  sequence  of  length  N  in  a  region  R 
containing  {  \. ( (A+B)_1A)}   if 


N-l 

Max 
zeR 

n    (1-t  z) 

n=  0 

=    min 
lTnJn=0 

Max 
zeR 

N-l 

n   (i-x  z) 

n       n 

n=  0 


More  conveniently,  we  will  say  the  best  sequence  in  R. 


In  the  remainder  of  this  section  a  best  sequence  of  parameters 
will  be  derived  for  a  circular  region  and  then  for  an  elliptic  region. 


22 


To  begin  the  derivation  the  following  definitions  are  required. 


Definition  3.2:   Let  g  (z)  be  any  complex  polynomial  of  degree  n,  then 
M  \z   (r  )1  denotes  the  maximum  modulus  function  of  g  (z)  for  |z|g .  r  ,   i.e., 


M  [g  (r)]=  Max 
n     izU  r 


gn(z) 


Definition  3.  3:     S     is  the  set  of  all  polynomials  g     of  degree  n  such 

that  g    (1)    =  1. 
ton 


Definition  3. k:     M  [g  (r,d)]  = 


Max 
d-  z  I  g  r 


gn(z) 


Definition  3. 3:  S  is  the  set  of  all  polynomials  g  of  degree  n  such 

that  g  (0)  =  1. 
n 


* 


The  polynomials  which  are  required  in  (3-1^)  belong  to  S   .  However 


n 


we  will  need  the  following  theorem,  due  to  E.H.  Zarantonello  (see  Varga 


[19]),  which  is  stated  for  polynomials  in  S  .   Therefore,  after  stating 
the  Theorem  it  will  be  modified  in  Theorem  3.2  to  apply  to  polynomials 

in  S  . 

n 


Theorem  3.1:     For  all  r   such  that  0  g  r  ^  1   ,     min  {M[g    (r)]] 

g  eS     n 
ton  n 


=  r  for  all 


positive  integers  n. 


23 

Theorem  3.2:   For  all  real  r  and  complex  d  such  that  0  ^  r  <  |d|, 


mi 


g  eS 
&n  n 


»  (M[gn(r,d)],  «  -£, 


for  all  positive  integers  n. 


Proof:  Let  r  and  d  be  such  that  0  ^r  <  |d|  and  suppose  there  exists 


some  g  eS  such  that  M  fg  (r,d)l< 
n  n  L  n     J 


Let  h  (z)=g  (d-dz).   Then 
n  '   n 


r 
T 


n 
> 


Max 
d-  z  I  ^  r 


gn(z) 


Max 
d-  z  I  ^  I  r 


M) 


z  g 


Max   |h  (z)  I   =  M[h  (  -£■  )]. 
1  n  n   d 

r 


Since  h  (l)  =  g  (o)  =  1,  the  above  inequality  contradicts  Theorem  3.1, 


and  M  [g  (r,d)]sL«j- 


for  all  g  is  S  .   On  the  other  hand,  if 
n     n 


Sn(z)  =("^r)n  thenM^[gn(r,d)]  = 


d-z  \  n 

theorem  follows 


and  the  conclusion  of  the 


n 


Corollary  3-2.1:   The  polynomial  P  (z)  =  ( — ■— — j   is  the  polynomial  in 

"X-  -&  -#- 

S  satisfying  M  [P  (r,d)]  =  min„  M  [g  (r,d)]  for  all  integers  n>  0  and 
n  n  _"*"      n 

g  eS 
n  n 

all  r  and  d  such  that  0  ^  r  ^  I d I . 


Proof :   This  follows  immediately  from  the  proof  of  the  theorem. 


Corollary  3.  2.1  can  be  used  to  determine  a  best  sequence  of  t  's 


to  use  on  a  disc  D  =  {z:|z-d|  <  r,  where  0  <  r  <  |d|]  .   That  the 


2k 

\.((A+B)~  A  lie  in  such  a  disc  is  seen  as  follows 


1 


The  idea  of  choosing  a  sequence  of  t  ' s  is  to  speed  convergence 
of  iteration  (2.5).  It  is  reasonable  to  assume  there  is  a  constant  value 
t  for  which  the  iteration  is  convergent.  From  this  and  (3-3)  it  is  seen 
that  for  some  t  >  0, 

X± ((A+B)"^)  e  B(-~  ,05-i-)-  {z:  ~~  -1  *  -i-}  -  (3-15) 

Consequently  Corollary  3.2.1  can  he  used  to  determine  a  best  sequence  of 


t  ' s  in  a  disc  as  follows, 
n 


Theorem  3.3:   If  the  eigenvalues  of  (A+B)"  A  satisfy  (3*15) 

then  the  best  sequence  of  t  's  to  use  in  (2.5)  is  t  =  x  a  constant. 

Moreover  t  =  — —  where  d  is  the  center  of  the  smallest  disc  containing 

all  X.  ((A+B)"^)  . 


Proof:   Since  no  confusion  will  arise,  X.   will  he  used  to  denote 
1 

X.((A+B)~  A)   in  this  proof.   Let  s(d,0;r)  he  the  smallest  disc  containing 

all  the  X. .   That  the  smallest  disc  is  of  this  form  is  seen  as  follows. 

1 

Since  the  complex  eigenvalues  occur  in  complex  conjugate  pairs,  the 
smallest  disc  containing  them  has  a  real  center,  (d,0).   In  addition 
(3.15)  implies  0  g  r  <  d.   Now  from  Theorem  3.2  and  Corollary  3. 2.1  it 
is  seen  that  of  all  polynomials  in  S   ,  the  minimax  polynomial  on  S(d,0;rj 

13  .H 


v*>  =  -¥ 


and  M*[PN(r,d)]  = 


N 

r 
d 


It  is  apparent  that  either  at  least  k   of  the  X.    lie  on  the 


25 
circumference  as  in  Figure  3- la,  at  least  2  as  in  Figure  3.1b  or  3  with 
one  on  the  real  axis  and  two  as  in  Figure  2.1a. 


Figure  3.1 


Let  s(&+€,  ,6pjr')  "be  any  other  disc  containing  all  the  X..      The 
polynomial  g^S  which  is  the  minimax  polynomial  on  this  disc  satisfies 


M*[qN(r',(d+ei;e2))] 


[(a^r+Cegn  "2-j 


N 


(3-17) 


The  Theorem  will  be  established  if  it  is  shown  that  the  equantity  in  (3.I7) 
is  greater  than  the  quantity  in  (3.16).  To  establish  this  we  shall  assume, 
without  loss  of  generality,  that  e,  and  ep  are  non-negative.   If  this  were 

not  the  case,  the  following  argument  would  be  changed  only  in  the  selection 
of  the  X.    of  Figure  3.1. 

Since  5(d+e ,ep;r ' )  contains  all  the  X.    it  in  particular  contains 
one  which  is  as  X   of  Figure  3.1a  or  Figure  3. lb, i.e.;,  X     =  (u,-v)  where 


u  >  0  and  v  ^  0.   Therefore, 
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(r,)2^(d.+e1-u)2+(e2+v) 


=  (d+€1)2+u2+e22+v2+2e2v-2u(d+€1) 

(d+€1)2+u2+€22+v2-2u[(d+e1)2+e22]~2" 
=   [((d+ei)2+e22)"2".u]2+v2. 
Hence  the  quantity  in  (3. IT)   is  greater  than 


{[(d+£l)2+e22]  2  -u}2+v2 
2" 


(**!>' 


+  e. 


N 
2 


Thus 

M*[qN(r',(d+e1+G2))]    > 

2       2 
where  e_=[(d+e_  )  +e0    ]  2     -d  >  0. 

3  J-         <= 


/*  \2     2 

(d+e   -u)   +v 

3 


(a+e3)< 


And 


M*[qN(rS(d+e1,e2))]   > 


(d-u)2+v2 


N 


=  M   [PN(d,r)]. 


The  last  inequality  is  derived  from  d  >  u  and  the  Law  of  Sines. 

We  have  proved  that  the  polynomial  PM(z)   of   (3.16)  has  the  minimum 

maximum  on  the  smallest  disc  containing     each       "K-  >  where  the  minimum  is 

taken  over  all  polynomials  in  S     ;    and  this  maximum  is  less  than  the  minimum 

maximum  of  all  polynomials  in  S     on  any  other  disc  containing  the  X..      The 

t     should,   therefore  be  chosen  so  that 
n 


N-l  .  .N  N 


i.e.;  t  =  — —  .      This  completes  the  proof. 
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If  more  is  known  about  the  X. ((A+B)"  A)  then  some  more  useful 
results  concerning  a  best  sequence  of  parameters  can  be  obtained.   In 
particular,  if  the  \. ((A+B)"  A)  are  contained  in  a  ellipse  then  the  best 
sequence  is  a  Tchebychev  sequence  as  in  the  real  case.   Before  proving 
this  we  need  to  recall  some  terminology 


2  2 

Definition  3.6:   In  the  ellipse   ^X"^   +  — «Z|J —  =  1  with  a2  s  b2, 

a  fe 

2   JL_ 
a  is  called  the  semi -major  axis ,  e  =  (l-  — -  )   2  is  the  eccentricity 

b 

and  (p,q)  is  the  center. 


A  modification  of  the  following  result  due  to  Clayton,  see 
Wrigley  [23],  will  be  used  in  determining  a  best  sequence  of  x  's 
when  the  X.((A+B)  A)  are  contained  in  an  ellipse. 

Theorem  3-k:      The  polynomial  in  S„  which,  over  the  ellipse  of  semi -major 
axis  a  <  1,  eccentricity  e  and  center  at  the  origin,  has  the  smallest 
maximum  modulous  value  is  given  by 

PN(z)  =   TN(z/ae)/TN(l/ae)  .  (3-18) 

If  the  center  of  the  ellipse,  still  assumed  to  lie  within  the  unit  circle, 
is  translated  to  (d,0)  then  the  minimax  polynomial  is 

PN(z)  =  TN((z-d)/ae)/TN((l-d)/ae).  (3-19) 

Remark :   As  e  ->  0  the  ellipse  becomes  a  circle  and,   as  expected,  the 
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polynomial  P^  of  (3. 18)  approaches  zN  .   This  is  easily  seen,  for  when 
e  is  small , 

T  (z/ae)  =  -|-  [exp(Ncosh"1(z/ae))+exp(-Ncosh'  (z/ae))] 


=  -i-  exp(Ncosh"  (z/ae)) 


1  -I  N 


=  4"  [IT  +(Wae)2-D-] 


.  2N-1(z/ae)N  . 

Hence  lim  P(z)  -   ^"fo"^     ■  zN  .   In  the  first  equality  the  fact 
e-*Q  iN       2   (l/ae) 

that  T  (z)  -  cosh(Ncosh'1z)  is  valid  throughout  the  complex  plane  was  used 
(see  for  example  Wrigley  [23]  ). 


* 


Since  we  are  interested  in  polynomials  in  S^  rather  than  S^  , 
Theorem  J>,k   is  rephrased  as  follows. 


Theorem  3.5:   Let  E  be  the  ellipse  with  semi -major  axis  a,  eccentricity 

N 


e  and  center  (d,0)  where  a  <  d.   The  polynomial  in  S  which  has  the 


smallest  maximum  modulous  value  on  E  is  given  by 

m   I       d-Z   \ 

TN(-^~ ) 

PN(z)  =  T  (d/ae)    '  (3-20) 

N 


Proof:   First  note  that  P  (z)  of  (3-20)  is  in  S*  .   Now  suppose  there  is 
some  qN(z)  in  S  with  smaller  maximum  modulous  than  P^  on  E.   Then 


Max 
zeE' 


qN(d-dz) 


Max 
<       zeE' 
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PN(d-dz) 


(3.21) 


where  E*  is  the  ellipse  with  semi -major  axis  a/d  ,  eccentricity  e  and 
center  (0,0).   But 

pH(d.dz) .   Jalzkil  ■ 

e  d 

Thus  (3.21)  contradicts  Theorem  3-h   and  the  conclusion  of  the  theorem 
follows . 

The  following  corollary  shows  that  polynomial  PN(z)  in  (3-20) 

can  be  put  in  the  form  of  the  polynomial  (3. 12).   Thus  the  best  sequence 
of  t  's  in  an  ellipse  can  be  implemented  in  the  same  way  that  the  best 

sequence  on  an  interval  (e  ,e  )  is  implemented. 


Corollary  3. 5. 1:   Let  E  be  as  in  Theorem  3» 5  and  suppose  each 
\j((A+B)  A)  is  in  E.   The  best  sequence  of  parameters  to  use  in  (2.5) 
on  E  is  the  sequence  {x  }  for  which 

N-l  TN((e1+e2-2z)/(e2-e1)) 

ni0   d-V}  =PN(2>  =   TN((e1+e2)/(e2-ei)) (3-22; 

where  e   =  d-ae  and  e   =  d+ae. 


Proof:   Let  E  be  as  in  Theorem  3-5-   As  noted  when  considering  a  disc, 
the  existence  of  such  an  E  which  contains  each    X.    is  guaranteed  by 
(3.I5).   Let  e   =  d-ae  and  e   =  d+ae.   The  best  sequence  of  parameters  on 
E  yeilds  a  polynomial 
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V>-JSo  (1-v)  =rrx-   *  (3.23) 

IT  ae  ; 
But  e  +e  =  2d  and  e  -e   =  2ae  so  that  (3.23)  is 


g^V^f   )     /e2+er2z 


N  '        /  Vel  \        .   /  e2+el 


e2"ei  J  N  V  62  "ei 


as  claimed. 


Remark :  As  the  ellipse  degnerates  to  a  line,  e  -*  1  so  that 


e  -*  d-a  and  e  -*  d+a 
which  is  exactly  the  result  of  the  previous  section. 


k.      USING  A  SEQUENCE  OF  t's  WHEN  LITTLE  IS  KNOWN 
ABOUT  THE  EIGENVALUES  OF  (A+B)_1A 

lf.1  Introduction 

In  this  chapter  we  will  develop  an  always  convergent  algorithm 
to  solve  Ax=q.   The  algorithm  is  similar  to  the  method  presented  in  section 
3- 3  hut  is  superior  since  it  does  not  require  any  knowledge  of  the 
eigenvalues  of  (A+B)_1A.   The  theoretical  basis  is  derived  in  sections  k.2, 
k.J>   and  k.k.      The  algorithm  is  stated  in  4.5  and  some  experimental  results 
are  given  in  4.6. 

The  sequence  of  t  's  in  the  iteration 

(A+B)xn+1  =  (A+B)xn  -  Tn(Axn-q)  {kml) 

is  defined  during  the  execution  of  the  algorithm.  As  in  the  methods 
presented  in  Chapter  3,  the  sequence  is  dependent  upon  an  interval.   Instead 
of  using  a  constant  interval  as  the  previous  algorithms  do;  the  algorithm  of 
this  chapter  generates  a  sequence  of  intervals  (a.,  b.)  which  converge  to  the 
optimal  interval. 

The  following  is  a  brief  description  of  the  algorithm.   It  is 
shown  that  5n  =  VXn-l'  which  is  automatically  calculated  in  the  practical 
implementation  of  {k.l) ,   can  he  made  to  approach  an  eigenvector  of  (A+B)_1A 
corresponding  to  an  extreme  eigenvalue.   This  approximate  eigenvector  is 
used  to  compute  an  approximation  to  the  corresponding  eigenvalue.   The 
interval  (a^  Td±)   is  then  improved  by  replacing  one  endpoint  with  the 
approximate  eigenvalue. 
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The  algorithm  to  be  developed  can  he  applied  to  any  iteration 
of  the  form  (l+.l)  in  which  A  and  A+B  are  positive  definite.   If  A  is 
derived  from  the  discretization  of  an  elliptic  partial  differential 
equation  as  described  in  Chapter  2  and  A+B  is  defined  by  one  of  the 
factorizations  due  to  Kim  [8],  Stone  p_5  ],  or  Dupont,  Kendall  and  Rachford 
[5]  then  A  and  A+B  will  be  positive  definite.   Thus  the  algorithm  can  be 
applied  to  these  factorizations. 

It  is  assumed  in  the  sequel  that  A  and  A+B  are  positive  definite. 
Additional  Notation;  Throughout  this  chapter  the  subscripts  and  super- 
scripts will  denote  the  following: 

i  corresponds  to  the  interval  (a.,  b.). 

n  and  N         correspond  to  the  number  of  iterations. 

k  and  K  correspond  to  the  dimension  of  the  matrix  A. 

Thus  the  vector  v  calculated  on  the  n   iteration  using  the  interval  (a.,  b.) 

will  be  denoted  by  v  .   This  vector  will  be  denoted  by  v  when  the  interval 
J     n  n 

used  is  not  being  stressed. 

The  algorithm  will  use  vectors  which  are  calculated  during 
execution  of  computer  programs  written  to  solve  the  set  of  equations  Ax=q 
by  using  a  factorization  and  iteration  (4.1).   Therefore,  we  shall  begin 
with  a  description  of  the  practical  utilization  of  the  methods  presented  in 
Chapter  3. 

Two  methods  are  presented  in  Section  3.3.   In  the  first  a  vector 
Xq,  an  interval  (a,b)  and  an  integer  N  are  selected.  The  iterates 
x  .  n-1,  . . .  ,   N,  are  then  defined  by 

(A+B)xn+1  =  (A+B)xn  -  -rn(Axn-q)  (4.1) 
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where  the  T  are  selected  so  that 


pw(x)  =      (i-T  x)  =  Nl  b-»     ;   . 

N      nio     n     T  (£*) 
In  the  second  method  the  iterates  are  defined  by 


x  , ,  =  x  +5x  . 
n+1    n   n 


5x   „  )Tn(y)     ,(A+B)-1(q-^  )  +   ^44s^  „  (4.2) 

85o=Jf(a+b)"VaS0), 

where  y  =  4; r  and  x_  is  the  initial  vector. 

(b-aj      0 

Throughout  the  remainder  of  this  chapter  x  will  denote  a  vector 
defined  by  (4.1)  and  x  will  denote  a  vector  defined  by  (k.2). 

The  vectors  in  the  following  definitions  are  calculated  in  the 

practical  use  of  (4.1). 

Definition  4.1:  The  n   residual  is  the  vector  r  =  q-Ax  . 
n   ^   n 

Definition  4.2:  The  n incremental  vector  is  8  =  x  -x  , . 

Using  these  definitions  the  iteration 

(A+B)x  ,_  =  (A+B)x  -  t  (Ax  -q) 
x   J   n+1   v   '   n    n   n 

can  be  rewritten  as 

(A+B)6  =t   ,r  .. 
'  n    n-1  n-1 

This  equation  is  solved  by  replacing  A+B  by  LU  and  calculating  5  from 

Lz  =  t   ,r  . 
n-1  n-1 

and 

U5   =  z. 
n 

The  next  iterate  is  then  x  =  x  -,  +  P> 

n    n-1   un 


The  vectors  x  are  similarly  calculated.   Instead  of  6  and  r 
n  n     n 

however,,  we  define 

A       „A 

Definition  k.J>:     ?n   =  q  -  ^n- 

Iteration  (U-. 2)  is  performed  hy  computing 

-1a 
z  =  L  r   , 

n  ' 

*n+l  -  U"lz  > 
and  ,A      ^Tn(y)    a     Vl^)  a  ^.5) 

6xn  =  (b-a)Tn+1(y)  Vl  +  T^yJ  5Vl 

n  >  1. 

The  vector  Sx»  is  computed  from 

t-1a 
z  =L  rQ 


and 


5xn  =  — -=-  U  z  . 
0   a+b 


For  every  n  >  0   x  ,  is  calculated  from 


A        A       A 

x  , ,  =  x  +  8x 

n+1    n     n 


Note  that  6  defined  in  (h-.~*>)    satisfies 
n 

,-1a 


The  importance  of  this  vector  will  he  shown  in  Section  4.4. 


4.2  Calculating  An  Eigenvector  And  Eigenvalue  Of  (A+B)~  A 

In  this  section  it  is  proved  that  if  y0  is  any  vector  and  if 

%  =PH((A+B)-1A)y0 

Wh6re  T  I  a+b-2X 


v«  -,*:;■' 


N  Ib-a 
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and  a>  X    .   or  b  <  \ ..._,  then  y.T  approaches  an  eigenvector  of  (A+B)~  A. 

mm        MAX       N 


It  will  also  be  proved  that 


<AyN,yN> 


WN       <(A+B)yN,y  > 
approaches  the  corresponding  eigenvalue. 


(4.6) 


-1, 


In  the  next  section  6„  ,  is  shown  to  satisfy  5  .  =  P  ((A+B)~  A)S,. 
From  this  it  follows  that  8  .  is  approaching  an  eigenvector  of  (A+B)~  as  N 


increases. 


We  shall  begin  with  some  motivation. 


Suppose  that  yn  is  any  initial  vector  and  that 


N-l 


'N 


(l-x*(A+B)"1A)y0  =  P*((A+B)_1A)y0,  where  P*(x)  =  II  (1-t*X).  (k.l) 


n=0 


N-l 

n 

n=0 


The  object  is  to  select  the  x     such  that  y  approaches  one  of  the  extreme 

eigenvalues  of  (A+B)~  A.   If 

K 
y  =  Y    avwv  *  where  the  w  ,   k=l,  ...,   K,  form  a  complete  set 
k=l 

of  eigenvectors  of  (A+B)"  A,  then  y  is  given  by 


K 


'N 


L  °k] 


k=l 


N 


^((A+B)"^) 


w. 


=  P. 


N 


^((A+B)"^) 


I    a      ^[^((A^B)-^)] 


k=l 


To  make  y  most  rapidly  approach  w  ,  {t  )in  (k-7)   should  be  selected  so  that 


Max 


P;[xk((A+B)-1AJ 
P;[X1((A+B)"1A)] 
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is  minimized.   Since  \0,    . ..,  \K  are  unknown,  we  minimize 


Max 

\K<  x<*.2 


*J« 


F*  [V^Br'A)] 


(^.8) 


where  \,  denotes 


\k((A+B)_1A). 


-1, 


Note  that  y  converges  to  an  eigenvector  of  (A+B)"  A  provided 


n 


lim        Max     — : — £-p r — r 

N^oo     k=2,...,K     PN  I  ^((A+B)"  A)J 


=  0  , 


a.9) 


or 


lim 

N-*oo 


Max 

\  <  x<  \2 


4<V 


=  0 


Minimizing  (1+.8)  merely  speeds  the  convergence. 

The  polynomial  P  ,  which  minimizes  (lj..8),  is  a  scaled  and  translated 
Tchebychev  polynomial  as  seen  by  the  corollary  to  the  following  theorem. 


Theorem  h..  1 :  Let  M(P)  =  Max 

-1<X<1 

degree  N  has  the  property  that 


P(X)|.   The  Tchebychev  polynomial  T  (X)  of 


for  all  N  degree  polynomials  q  and  all  y  jt   [-1,1] 


Proof :  The  theorem  is  a  restatement  of  the  classical  result  of  Tchebychev 


that  M(T  )<M(q  )  for  all  q  such  that  CLT(l)=CL,  and  |T(y)  >  q__(y) 


*B' 


m 


IV 


'T 


*W 


for  all 


y  /[-1,1]. 
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Corollary  U.l-1:   The  polynomial 


P*(X)  = 


I   a+b-2X  ] 
Nl  b+a   / 

N  I  b-a  I 


satisfies 


Max 
a<  X<b 


P*(X) 


Max 
a<  X  <  b 


qN(X) 


PN(y) 


< 


qN(y) 


for  each  N  degree  polynomial  q  such  that  q  (0)=1  and  for  all  y4  [a,b]. 


Proof:   Omitted. 


Theorem  1+.2;   Let  w  ,  k=l,  .  . .  ,  K,  form  a  complete  set  of  eigenvectors  of 
(A+B)~  A  corresponding  to  the  positive  eigenvalues 


\.,fl v  =  X     >  X0  >  .  . .  >  Xv  i    >  K  =  x    •   >  °- 
MAX    1    2  —  —     K-l    K    min 


K 


Let  yn  =  V  a,w  be  any  vector  such  that  c*  a  /:  0  and  define  the  sequence 
k=l 


of  t    .   n=0,    ....   N-l.    such  that 


a+b-2X 


N-l  b-a 

p  (x)  =  n  d-xnx)  = 


n=0 


t(— ] 

Nlb-a/ 


(U-.10) 


where  a>  X    .      or  b  <  X...„  and 

min  MAX 

then  either 


Nv   mm7 


^ 


VHlAX^ 


If  %=V(A+B)"lA)yo 


lim  1 


Nv   MAX' 


yN  =  Vl 


or 


(4.11) 


S!     r   ;'       .      yN=VK 


38 


The  convergence  of  yN  to  w1  or  wK  is  most  rapid  when  a  =  \^n 
and  b  =  L  or  when  a  =  X^  and  b  =  X^,  respectively.  And  y^   converges 
to  w  only  if  b<  Km^   and  converges  to  wK  only  if  a  >  X^^. 


Proof:  Suppose  a>  X    .   or  b  <  \ 
,t  **  mm 


MAX' 


Suppose 


vw 


< 


"VSlAr 


This  implies  b<  \,,.^  and 
MAX 


'H 


(XMAX} 


> 


?N(V 


k  =  1,  . .. ,  K. 


From  the  definition  of  yQ  and  y^  we  have 

i  *      w 

— t  %  =  L. " 


^T^MAX' 


i"k     VW        "k 


a+b 


K 


-»* 


■E 


a, 


"Nl       b-a         j 


k=l         T. 


N 


a+b -2 


:>*MAX 


=  aiWl  +  .     J 


b-a 


a. 


a+b -2^ 


•N 


b-a 


w. 


vra,»i*TJa+b-aWl    " 

I        b-a 


a+b-2\.  1 
k 


)  TNl      b-a      j 

St?    .-,     k     m   /a+b-2\ 


>^[a,b]  ^  TN 


SlAx] 


w. 


b-a 


(4.12) 


Applying  Corollary  if.  1.1  to  the  terms  of  the  first  sum  we  see  that  as  N 
increases  these  terms  are  decreasing  at  an  optimal  rate. 

For  the  terms  of  the  second  sum  we  have  from  (if.  12)  that 


1  < 


a+b-2\. 


k 


|b-a| 


< 


a+b-2\. 


MAX 


b-al 


59 


Thus 


cosh 


-1 


a+b-2\, 

I 

b-a 


-   cosh 


-1 


a+b-2\. 

b-a 


MAX 


<  0. 


Using  the  inequalities   and  the  relation 

T   (X)    =   [exp(N  cosh_1X)  +  exp(-N  cosh_1X) ]/2 


W 


for     X     >  1     yields 


lim 


N 


a+b-2\ 


k 


b-a 


=  0 


Tja+b-^MAx' 


b-a 


Thus  lim  Tr-77 r     y*T  =  a-1wn    which  establishes  the  first  half  of   (k.ll). 

The   second  half  of   (k.ll)   is   similarly  derived  assuming     P,T(\    .    )   >  P,T(X„„„) 

N  min  J    N  MAX 

That  the  convergence  is  most  rapid  when  a  =  \   .      and  b  =  \_  or 

min         2 

a  =  \  ,  and  b  =  ^MAY  follows  directly  from  Corollary  4.1.1. 


Finally  we  note  that  y, ,  yp,  ...  converges  to  w  only  if 


V^mat 


> 


V\iin} 


to  w_.  only  if  a  >  X   . 
K    J        mm 


,  i.e.,  only  if  b  <  X 

This  completes  the  proof. 


MAX.   Also  y1,  j2,    ...  converges 


Any  approximation  to  an  eigenvector  of  (A+B)"  A  yields  a 
corresponding  eigenvalue  approximation.   This  is  derived  from  a  theorem  in 
[20],  presented  here  for  convenience. 


Theorem  L.  5  [20]:   If  A  and  B  are  positive  definite  then  the  eigenvectors, 
x,  ,  and  the  corresponding  eigenvalues,  \,f   of  the  generalized  eigenvalue 


problem, 

satisfy  the  following: 

1,  \  >  0; 

2,  xkV.  =0   j  /k. 


A^  =  XkBXk  , 


ko 

Proof:  Since  B  is  positive  definite  there  exists  an  orthogonal  R  such  that 

R  BR  =  diag(&  )  =  D  ,   where  Pk  is  an  eigenvalue  of  B. 

(A-XB)  =  RD(D"1RTARD"1-XI)DRT. 

Hence  det(A-XB)  =  (det  R)2(det  D)  det(P-Xl),  where  P  =  D-1R  ARD-  .   Since 
(det  R)  =1  and  (det  D)  =  H(p,  )  >  0,  the  zeros  of  det(A-XB)  are  those  of 
det(P-Xl).   These  are  the  eigenvalues  of  a  symmetric  matrix  and  are  therefore 
real.   The  matrix  P  has  a  complete  set  of  eigenvectors  z,  which  may  be  taken 
to  be  orthogonal.  Hence 

D_1RTARD_1z  =\  z  ;  and  A(RD_1zk)=?^RDzk  =  \kRD(DRTRD"1)zk  = 
X,  B(BD"  z.  ).   Thus  x.  =  RD~  z,  is  an  eigenvector  of  the  generalized  eigenvalue 
problem,  corresponding  to  X,. 

Since  {z  }  is  orthogonal  and  RD"  is  nonsingular,  {x,  }  forms  a 
complete  set  of  real  eigenvectors.   Furthermore 

m  m     m    rn  mm  m 

0  =  zk  z.  =  (DEx,)  DBx.  =xk  RBDRx.  =  x^  Bx .  , 

so  that  {x  }  is  orthogonal  with  respect  to  B. 

Finally,  X,    >   0  since  <  Ax,  ,   x,  >  =  <\,  Bx,  ,  x>   This  completes 
the  proof. 


Observe  that  this  theorem  applied  to  A  and  A+B  implies  that  the 
eigenvectors ,   w ,   of  (A+B)~  A  are  orthogonal  with  respect  to  (A+B). 
The  next  lemma  defines  an  approximate  eigenvalue. 


K 
Lemma  l+.l:  Suppose  A+B  and  A  are  positive  definite.   Let  y  =  *£,     e\r\ 

k=d 

where  wk  is  an  eigenvector  of  (A+B)~  A  corresponding  to  X ,  .   Let 

U   =    NYjAy  >    .   If  y  is  approximately  equal  to  w  with  errors  €  in  the 
<y,(A+B)y>  J  K 


kl 


direction  w ,  ,   and     e     »  e 

^2 


k  j=  J  _,  then  \i  =  X       with  error  of  order 

J 


Proof:  Assume  J  =  1.   The  remaining  cases  are  similar.  Assume  {w  }  is 
orthonormal  with  respect  to  (A+B).   Let 


M  ~  <y,  (A+B)y> 


U.13) 


K 


Substituting   ]T  e  w  for  y  and  (A+B)(A+B)"  A  for  A  in  (4.13)  yields 


k=l 


"k  k 


W  = 


K 


K 


2_  €vwvT(A+B)(A+B)'1A  ;  ,  e,w. 


k=i 


"k  k 


I 

J-l 


J  J 


K  K 

k=d  j=l 


€  .W  . 


1  T 

Now  using  the  facts  that  (A+B)"  Aw.  =  \.w.  and  w  (A+B)w-  =  5    ,  where 

J  J    J  .k  J  j  ^k 

5.  .    is  Kronecker  delta,  we  obtain 

Al  k    k  1     k^2  k     k     X 

£       2  "  £        2/   2 

Z  \       1+L  v€i 

k=l     k  k=2    k     X 


This  completes  the  proof. 


In  the  case  A+B=I  the  above  lemma  reduces  to  a  result  of  Hestenes 


[21]. 


k*3     Application  Of  The  Method 


To  compute  an  eigenvalue  using  the  above  method  it  is  necessary 
to  determine  a  sequence  of  vectors  y  ,    ...,   y  such  that  y  =  P((A+B)~  A)y 


where 


VX)  =  TH 


a+b-2X 

b-a 


"N 


a+b 


b-a 


k2 

The  next  theorem  shows  that  8  ,   n.  =  1>  . ..,  N  +  l  is  the  desired  sequence 

of  vectors.   Since  s  is  calculated  when  iteration  (4.1)  is  used  to  solve 

n 

Ax  =  q,  no  additional  work  is  required  to  obtain  the  approximate  eigenvector. 


Theorem  4.4;  Let  x0  be  any  vector  and  define  x  ,  n  =  1,  2,    ... ,   N+l  by- 
iteration  (4.1).   If  the  sequence  of  fr  )used  in  the  iteration  is  such  that 


NJ.  T  A^m) 

P^(X)  =  nQo  (1  "  ^  = 


N  I  b-a  i 


where  a  >  \    .      or  b  <  \„AV  and 

mm        MAX 


PN(Xmin) 


/ 


PN^MAX' 


'  then  Vi  =  Vrx 


N 


approaches  an  eigenvector  of  (A+B)~  A  as  N  increases. 


Proof:   Since  8n  =  xn  -  x^  =  (i  -  ^(M^A)*^  +  Vl^'W  > 


,-1 


-1, 


5  -t   ,(A+B)"q=-T  _(A+B)  Ax  . 
n  n-lv   '  H     n-lv   '    n-1 


(4.U) 


and 


5  ^.-t  (A+B)"1q  =  -t  (A+B)_1Ax  . 
n+l  nv   /  m,     n^   /    n 


(4.15) 


Subtracting  (t   /t  )  times  (4.15)  from  (4.14)  yields 


n-1'    n 


5-(x     ,/t    )8  ^n    =t     .  (A+B)'1A(x  -x     n)    =t      n(A+B)_1AS     , 
n       n-1'    n'  n+l         n-lv        '        v  n     n-1'         n-lv        '         n  ' 


k-l. 


which  implies   (I  -  x^A+B)"  A)8n  =  (Tn_l/Tn^5n+1  * 

Hence    (Tn/Tn-l)(l"Vl(A+B)"lA)5n   =  5n+l'      And 

N-1 
VT0       IL   (I   "   Tn(A+B)"  A>51   =  Vl     ' 


Since  the  sequence  of  t  ,  n  =0,  .  . .  ,  N-1,  is  used  cyclically,  T  =  t 


+3 

and  thus 

N-l 
5N+1  =      (I-t  (A+B)-  A)5l  .  (4.l6) 

n=0 

The  conclusion  follows  from  Theorem  4-2. 

An  approximation  to  the  corresponding  eigenvalue  can  be  calculated, 
according  to  Lemma  4.1,  by 


Mm  "  <WAtB%i> '  (     ' 


4.4  Using  The  Recursive  Property  Of  The  Tchebychev  Polynomials 

From  the  theory  developed  in  Sections  4-2  and  4.3  it  is  possible 
to  use  iteration  (4.1)  to  solve  Ax=q  and  simultaneously  compute  an 
approximate  eigenvector  of  (A+B)~  A.   The  process  is  started  by  selecting 
N  and  using  (4.1)  to  compute  the  sequence  of  vectors  x  ,  n=l,  . .  . ,  N+l. 
The  resulting  approximate  eigenvector  is  5^ _  =  x  , -x_,  .   In  the  next 
theorem  it  is  shown  that  if  x  is  computed  from  (4.2)  then 

%n+1    =  (A+B)-1(q-A^n)  (4.18) 

is  an  approximate  eigenvector  for  every  n  =1,  2,    ...  .   Recall  that 

§  ,,  /x  , ,  -x  .  Note  that  §  n  is  computed  when  (4.2)  is  used  to  solve 
n+l    n+l  n  n+l 

Ax=q  and  thus  the  approximate  eigenvector  is  again  obtained  with  no  additional 
work. 


Theorem  4.  5:  Let  x  be  any  vector  and  define  the  sequence  of  vectors 

x  ,  n  =1,  2,    ...,  by  iteration  (4.2).   Assume  the  interval  (a,b)  used  in  the 


kk 

iteration  satisfies  (a/b)  i>  (X^,  >-mx)-      ^t   rn  =  q-Axn-   It  follows 
that  5  =  (A+B)"  r  approaches  an  eigenvector  of  (A+B)~  A  as  n  increases. 
Furthermore  5  approaches  the  eigenvector  corresponding  to  X^^  only  if 
b  <  X    and  approaches  the  eigenvector  corresponding  to  X  .   only  if 

JXlA-X.  iiij.ii 


a  >  X 


min 


Proof:     Define  t       n  =  0,    ...    N  - 1,   such  that 


a+b-2x 

iio       n         t.  r    a+b 


[     a+b     I        » 
L    b-a    J 


and  let  t     =  t   /      ,  __v  for  n  >  N.      Let  x_   =  xn  and  define 

n    n(mod  NJ      —         0    0 

W    N      W  N    N  /    \-l/       \ 

x  =  x  n  +  5    where   8  =  r     _  (A+B)   (q-Ax  ,), 
n    n-1    n  n    n-lv   '   VM-   n-1' 


The  vectors  £   defined  be  (4.2)  have  the  property  that 
n 


a     n 

x  =  x   for  every  n. 

n    n         ° 


Therefore  ^  =  (A+B)_1(q-Axn)  =  (A+B)_1(q-AxJJ)  =  -^(A+B)_1(q-Ax^). 


Ttf-1 

But  xn  =  t^  and  Tn(A+B)_1(q-Axn)  =  5n  , .   Hence 
n    0      nv  '     VH   zr    n+1 


'n 


From  (ij..l6)  we  obtain 


n 


Jn+1 


n     un+l  ' 

:o 

n-1 


l.\JB 


Jn+1 


n"     IL(l-Td(A+B)     AK 


(4.19) 


3  =  0 


Since 


n       wl 


8?=   -4 


T*    (A+B)"1(q-Ax^) 


,-1 


=   (A+B)-±(q-Ax0)    =  5\    , 


*5 
(4.19)  can  be  rewritten  as 


%n+1=      H  (I-t^A+B)-1^  .  (1..20) 

j  =  0      J 


The  conclusion  follows  from  Theorem  4.2. 


The  approximate  eigenvalue  is  now  given  by 

We  now  show  that  no  matrix- vector  multiplications  are  necessary  to 
compute  the  quantities  (A+B)s  _  and  Ag  -,  which  are  needed  in  the 
calculation  of  ju  .    The  vector  (A+B)^     =  ^n   and  therefore  is  available. 
That  A5  -  can  be  calculated  without  a  matrix-vector  multiplication  is 
proved  in  the  next  theorem.   More  precisely,  it  is  shown  that  AS  ,  may  be 
calculated  from  known  quantities  with  only  two  vector  additions,  at  the 
expense  of  additional  storage  to  hold  (£  ,  -  £  ). 


Theorem  4.6:   If  £  is  the  n   iterate  defined  by  (4.2)  and  $  =  q-Ax  ,  then 
- —       n  n   ^    n' 

Tn+2(y)(b-a) 
AK+1   =  4Tn+1(y)      L^n+l"^n+2^qn+l^rn"rn+l^J   ' 

where  y  = (a+b)/(b-a)  as  in  (4.2)  and  q^  =  (Tn(y))/(Tn+2(y)). 

4T   (y) 

Proof:  From  (4.2)  we  have  £  ,_  =  x  ,n  +  %   ,. /  w, — r  +  q  ^-,  5x  ,   where 

v   '  n+2    n+1    n+1  T  pKYiVo-a.)        nn+l  n' 

y  =  (a+b)/(b-a),  q  .  =(T  (y))/(T. p(y))  and  &     is  given  in  (4.2).  Therefore 


■n+2 


+  T_n(y) 


^n+2  =  »  -^n+2  =  Vl  "  AVl  V"(y)(b-a)   "AvA 


14-6 

Replacing  q        A§n   by  <attfl(A£n+1  +  <1-<1-A&n)  and  collecting  terms  yields 

4T      (y) 
V2  -  Vl  =  Afl     Tn^(y)(b-a)     "  Vl(VVl] 

and 

T^>(y)(b-a)        r 

AVl=    T^73 [<WW"WVW 

which  completes  the  proof. 


We  shall  now  develop  a  method  of  using  iteration  (4.2)  to  solve 

Ax=q  and  automatically  improve  the  interval  (a,b).   The  method  will  use 

iteration  (4.2)  and  a  sequence  of  intervals  (a. ,b. )  that  converge  to  the 

optimal  interval  (\  .  ,  \   .S).      The  method  is  developed  from  Theorem  4.5 

and  Lemma  4.1.  According  to  Theorem  4.5,  as  n  increases  %     approaches  an 

eigenvector  of  (A+B)~  A  corresponding  to  an  extreme  eigenvalue  and  u. , 

defined  in  (4.2l)_,  is  an  approximation  to  the  corresponding  eigenvalue. 

The  intervals  are  improved  by  replacing  a.  with  \i.    if  it  is  an  approximation 

to  the  smallest  eigenvalue  and  replacing  b.  with  \x.    if  it  is  an  approximation 

to  the  largest  eigenvalue. 

Whether  u.    approximates  \   .      or  approximates  X,..,.  is  determined 
1  mm     **  MAX 

as  follows . 

Recall  that  §  approaches  an  eigenvector  corresponding  to  \   . 

only  if  a.  >  \   .   and  approaches  an  eigenvector  corresponding  to  X„„^  only 
1    mm  jr     -o     MAX 

if  b.  <  \,,AV.   Thus  ju.  can  approximate  \   .   only  if  a.  >  X  .   and  jti.  can 

1    MAX        1      r*  mm    J  1    man      1 

approximate  T^y  only  if  b.  <  Tw™-.   Furthermore  we  are  developing  a  method 

which  will  generate  a  sequence  of  intervals  (a.  ,b.  )  that  converge  to  (X   .  A„AV.). 

1'  1  mm'  MAX 

Therefore  we  must  be  able  to  improve  the  left  endpoint  of  some  intervals  and 

the  right  endpoint  of  others.   This  implies  that  the  condition  a.  >  X   .      and 

1    mm 


hi 
b.  <  \j,AV   must  be  satisfied  for  each  i.   If  this  condition  holds  it  can  be 
concluded  that  ju.  approximates  \   .      if  ju.  <  a.  and  approximates  X         if 


ju.  >  b.. 
l—i 


Note  that  the  condition  a.  >  X   .      and  b.  <  \....r   can  be  replaced 

l    mm      i    MAX         * 


by  a.  >  X  .   and  b.  <  \..AV,  i.e.,  (a.,b.)c:  (\  .   \„.„).      With  this  relaxed 
J   i  —  mm     i  —  MAX         i'  i     mm'  MAX 

condition,  the  endpoints  of  each  (a. ,  b. )  which  are  not  already  optimal  can 

be  improved. 

The  next  lemma  shows  that  if  (a. ,b.  )  cz  (\  X.„„)   and 

i'   l  mm'        MAX' 

(a.    , ,b.    , )   is  defined  by  replacing  either  a.    or  b.    with  a. ,   then 
l+l '   l+l  11  i' 

(a.  ,_  ,b.  ,  _  )  c  (\   .    ,  \.,AV).      Therefore  if   (a-,  ,bn  )  c  (\  .    ,  \,„AV)  we  may 
v  i+I '   i+I  mirr     MAX  v  1'  ly        v  mm'     MAX  ^ 

assume    (a.  ,b.  )  c  (\   .    ,   \.,A„)   for  each  i.      How  to   select    (an  ,b,  )  c   (\    .    .\rliV) 

v  i'  l      mm'  MAX  1'  1'    v  min'  MAX 

is  discussed  in  section  k. 5* 

Lemma  1+.2;  If  ju  =  ^vTlpS 5.  where  A  and  A+B  are  positive  definite,  then 

HW(a+b>~1a>>W(a+b>"1a>]   • 

Proof:   This  is  a  restatement  of  (3.6) 

Throughout  the  remainder  of  this  chapter  P  (X)  will  denote  the 

polynomial  .,  —„.        ,, 

T   (    a:Tb-2X   )/T   (^L.),    a<b   . 
n       b-a        "    n     b-a   " 

The  method  we  have   developed  so  far  may  be   summarized  as  follows: 

1.      Given  a  vector  &_   and  an  interval    (a. ,b. )   compute  x       n  =1,    ...    N 

0  1'    1  n 

and  fc1     n  =1,    ...    N  -1  by   (4.2). 


2.      Calculate 


<4-i,<i> 
1     <li/A+Bd> 


k8 
3.   If  n±  <  a±   define  a±+1  =  H±   and  b±+1  =  b± 
It  u.   >  bi  define  ai+1  =  a±  and  b±+1  =  fi± 
In  either  case  let  x^1   =  x^J  and  return  to  step  1. 

We  now  consider  the  case  a.  <  ju.  <  b..   If  the  approximate 
eigenvalue  is  repeatedly  in  the  interval  (a.,b.)  then  one  of  the  following 
conditions  holds:  Either 

1.   (a.,b.)  ^  (^  Xmx) 


or 


2.   P_(\  .  )  = 

IT  mm' 


P]TXMAX' 


Remark:  We  can  choose  (a,  ,b,  )  c  (\  .  _,  \MAV)  and  thus  assure 

(a.  ,b.  )  c:  (\  .  ,  \,„__).  However  we  will  present  a  heuristic  method  of 
v  xf   i     nmr  MAX 

selecting  the  initial  interval  that  does  not  guarantee  (a.,  ,b  )  c  (\  .  ,\   ) 

(see  Section  4-5).   For  this  reason  it  is  necessary  to  consider  case  1  above. 

In  either  case  the  vectors  x  will  converge  to  the  solution  of  Ax=q.   In 

the  first  case  the  new  method  reduces  to  the  method  of  Chapter  3«   The 

initial  interval  can  be  chosen  so  that  if  this  case  occurs,  convergence 

is  rapid  (see  section  k-5)-      In  the  second  case  convergence  holds  since 

0  <  P  (X   .  )  <  1  for  all  n  (see  Fig.  4.1).   Nevertheless,  there  is  a 
n  man  v      °  '  ' 

complication. 

We  shall  consider  two  subcases. 

2a.   a.  =  X   .  and   b.  =  X,,,,.  . 

l    min         l    MAX 


2b.   a .  £  X  .  and   b.  ^  \.,A__ 

l  '   mm         i  '   MAX 


Consider  2a.   If  a.  =  X   .   and  b.  =  \„.„  then  KT(\   .    )  = 
l    mm      l    MAX      Nv  nmr 


PN^XMAX^ 


Thus  if  the  interval  is  improved  so  that  (a.,b.)  =  (X   .    ,   \„AV),  then 

v  i'  i'   v  min'  MAX" 


PN(Xmin) 


V^MAX^ 


In  this  case  the  interval  is  nearly  optimal  and 
should  be  used  in  the  subsequent  iterations  to  compute  AjL-.,  ^+2*  *  *  * 


k-9 


In  case  2b  the  situation  may  be  as  pictured  below. 


p„(«) 


Convergence  of  fx_Tl  will  be  slow  since  the  interval  (a.  ,b.  )  is  a  crude 
approximation  to  (X   .  ,  X        ).      Furthermore  the  vectors  5   do  not  approach 
an  eigenvector  of  (A+B)"  A.   Instead  they  approach  a  linear  combination  of 

Note  that  merely  increasing 
the  degree  of  P,  i.e.,  performing  more  iterations,  will  not  guarantee 


the  eigenvectors  corresponding  to  X   .      and  Xr,.„ 
to  *•    ~o     mm      MAX 


convergence  of  {  §  )  to  an  eigenvector  of  (A+B)~  A.   It  follows  that  the 
interval  cannot  be  improved  using  the  method  described  above. 

We  shall  now  describe  a  technique  which  will  force  convergence  of 

Performing  one  iteration  of  (4.1)  with 


[&-)   when  Pjx   .  )  = 
1  nJ      I\P  mm' 


PN^MAX' 


t =  l/b.  will  diminish  the  components  in  the  error  vector  and  the 
approximate  eigenvector  corresponding  to  the  larger  eigenvalues.   Therefore 
it  will  increase  the  relative  size  of  the  components  corresponding  to  the 
smallest  eigenvalues. 

This  idea  can  be  stated  more  rigorously  as  follows.   If  N+l 
iterations  have  been  done  using  the  interval  (a. ,b. )  and  &   =P ((A+B)~  A)s 


where  P„(X  .  ) 
N  min 


PN^MAX^ 


,   then  one  iteration  with  x  =  l/b.  yields 


Va  =  W  ((A+B^A^ 


*  (         \ 
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where  P^X)  =  (l  -£)  PjCX).   Furthermore  P^O^)  # 

We  could  next  compute  r^+1  =  q-Ax^.^ 

V2  ■  (A+B)"\i  • 

i     <  V2,AV2> 

and  improve  the  interval  if  u^  /  (ai,b  ). 

However,  if  u  p  e  (a.  ,T>.  )  we  assume  that  case  2a  has  occurred 
and  thus  want  to  continue  iterating  using  (a. ,b.).   There  are  two  methods 
of  using  the  interval  for  the  next  W  steps.   The  first  is  to  define 
xj+1  =^+1  and  compute  x^+1  n=l,  ...  N  using  (ai+1,  b±+1)  =  (a^K).   The 
second  is  to  continue  using  (a.  ,b.)  and  thus  define  £LpJ  ^i\r-V  """  ^PN+l 

We  take  the  second  alternative.   Furthermore  the  calculation  of 
jli__,  0  and  the  corresponding  attempt  to  improve  the  interval  is  omitted. 
Instead  ^L.,  is  defined  using  T=l/b.  and  ^.p;  •  •  •  .>  ^pivr+p  a2,6  immediately- 
computed.   Then 

i        Nb2N+l,Ab  2N+1 

is  computed  and  used  to  improve  the  interval  if    ^otvu-i  e  (a-.»b.)  .   This 

procedure  is  chosen  for  the  following  reasons.  An  approximate  eigenvalue  u 

is  computed  only  once  every  N  iterations  instead  of  twice.   Next,  we  know 

that  case  2a  will  eventually  occur  in  every  application  of  this  method  and 

this  case  is  not  benefitted  by  the  extra  calculation  of  uL0.      Finally  case 

2h  happens  only  infrequently.   (This  can  be  established  by  fixing  X   .   and 

considering  the  range  of  values  of  \,,._r  that  cause  case  2b  to  occur.  ) 

MAX 
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A  technical  problem  now  remains.   The  vectors  x   ,  ...  x^ 

cannot  be  computed  by  simply  restarting  the  iteration  using  (4.2)  since 

SxV.  is  not  defined.   The  following  will  elucidate  the  problem.   Computing 

x  p  requires  £  ,  and  6X  .   First,,  5x    is  computed  and  then 

x1  „  =  &  ,  +  5x  , ,  .   If  iteration  (4.2)  is  used  to  calculate  x^",  ...  x^l 
n+2    n+1     n+1  1'     TT 

then  the  last  vector  5x  to  be  calculated  is  8$  , .  If  £  is  computed 
by  using  t  =  l/b.  in  iteration  (4.1)  then  it  is  necessary  to  also  define 
5x7,  before  resuming  iteration  (4.2). 

The  vector  S^L.  must  be  defined  so  that  the  iterates  &  ?,  ^vr,,^  •••> 
converge  to  x  and  so  that  oVr+2^  %4-v  '•'>    converge  to  an  eigenvector  of 
(A+B)~  A.   Convergence  of  [zl    }  is  established  by  proving,  in  the  following 

theorem,  that  the  error  vectors  e   =  x  -  £   n  >  N+3  satisfy  e  =P  ((A+B)~  A)e_, 

'  n      n    —  ^  J     n   n         0' 

where  again  P  (X)  =  (l  -y—)?     -,  (X).   Then  using  the  fact  that  the  same 

i 
polynomial 

PUA+B)-^)  =   n  (I-t  (A+B)"^) 
2=0      J 


appears  in 


and 


e   =  P  ((A+B)  1A)e. 
n    n  0 


U=V<A+B)  A)^i 


(compare  (3.8)  and  (4.20))  we  conclude 


U  =  p>+b>~1a>\  > 


and  thus  8  approaches  an  eigenvector  of  (A+B)  A. 


Theorem  4.8:   Let  xn  be  any  initial  vector  and  let  x  n=l,  . . .  ,   N  and 
5xn  n=0,  ...,   N-l  be  defined  by  (4.2).   Let  x^+1  =  ^N  +  £^+1   > 
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where 


and  let 


Vi  =  ^b)"Va\) 


e£N  =  (i-tCa+b)-1^^^  . 


If  x     ,  n  >  N+2  is  generated  by 

a  a   ,     A 

x     .    =  x  +5x     , 
n+1         n       n 

^n  -     (b-a)Tn(y)      (A.B)"^^)  +  ^-  ^   ,  ft.*) 

where     y  =  (a+b)/(b-a)    .     Then  •  a+-b_2x 


'n  1 "b-a 

C 

n 

for  n  >  N+2. 


en  =  (l-T(A+B)"1A)Pn((A+B)-1A)e0    ,   where  Pq(x)    =       n  lg  ^ 


Proof:     The  proof  will  he  by  induction  one  ,   j=l,   2,    ...    .      The  error 

vectors  e     =  x-x       n=l,    ....   N  satisfy     e     =  P  ((A+B)"  A)e^     where 
n  n  '  '  ^        n         n  0 


T 

p  (x)  =.  -S 


I  a+b-2X 
b-a 


T 


n  I   b-a 

Since  T  +1(X)    =  2XT   (X)  -T  (x)   for  all  n,    and     y  =  (a+b)/(b-a) 

it  follows  that 

Tn(y)U  T        (y) 

K  4-1  (X)    =   -       rp          /      \/v        \       p    (x)  +    P    (x)    +    m,      .     (P    (X)-P      .  (X)). 

n+1'                 T       (y)(b-a)       nv  nv    '       T       (y)    v  nv    y     n-lv    " 

ni  n+1                                     (If.  23) 


Consequently  the  polynomial  P   (X)    =  (1-tX)P     1(X)   satisfies 

*      .    .  Tn(y)^  *  *  T        (y) 

^1(y)(b-a)   WX>+WX>  +  T^ 


£*W    =  -  T°     (y)(b-a)   4i»Kkl«  +  T^UJ^tl^K^'  ^ 
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-1, 


Let  Vl   =  \  +   ^N+l   '   Where   Vl   =  T(A+B)~    (ci-A^f)   and  let 


-1. 


gftjj  -  (I-tCCa+B)-^)^^. 

Obviously  eN+1  =  P^+1((A+B)_1A)e0.  (4.25) 

Now  suppose  x   n  =  N+2,  . . .  ,   is  defined  by  (if.. 22).   We  have 

A      _     A         A 

eN+2  ~  X  "  XN+2  ~  X  "  ^+1  "  6XN+1 

h-T   (y)  T   (y) 


But  (q-AxN+1)  =A(x-xN+1)  =  A(l  -  T(A+B)_±A)eN  and 

=  (l-r(A+B)-1A)(eN_1-eN) 
Therefore  (4.26)  can  be  rewritten  as 


eN+2  =  (I-T(A+B)_1A) 


L 


WN(y)         x     TN1(y) 
eN"  (b-ajT^y)  (A+B)  AeN  +  T^T?)  <•*-*-!> 


With  (4.23)  and  e  =  P  ((A+B)  XA)e   this  becomes 

eN+2  =  (l-T(A+B)"1A)PN+1((A+B)"1A)e()  =  P*+2(  (A+B)_1A)e0  .      (l|.27) 

From  (I4.25)  and  (4.27)  e    =  P*   ((A+B)_1A)e  for  j=l  and  j  =2 .   The 
general  inductive  step  follows  trivially  from  (4.24)  and  5x =  e   .   -e 


4.5  The  Algorithm  And  Its  Convergence 

In  this  section  a  method  of  selecting  the  initial  interval  (a. ,b. ) 
will  be  developed.   A  precise  definition  of  the  algorithm  to  use  (4.2)  and 


simultaneously  improve  the  endpoints  a.  and  t>.  is  given.   Convergence  is 
also  shown.   In  the  next  section  some  empirical  results  -will  be  given. 

Effective  application  of  the  methods  of  the  previous  section 
requires  the  selection  of  an  initial  interval  (a^t^)  c  (^j^Wy). 
Lemma  4-3  (below)  implies  that  the  only  interval  guaranteed  to  satisfy 
this  condition  for  every  problem  is  the  point  interval  a,  =b.  =1.  Although 
the  methods  of  the  previous  section  work  when  this  point  interval  is  used 
as  the  initial  interval,  the  resulting  iteration  may  start  inefficiently. 
Indeed,  the  error  vectors  e  may  even  grow  until  the  interval  is  improved. 
Therefore  the  iterations  performed  after  the  interval  is  improved  must 
annihilate  an  error  vector  that  is  larger  than  the  initial  error  vector. 
We  shall  now  consider  an  alternate  method  of  selection  to  overcome  this 
inefficiency. 

Our  goal  is  to  develop  a  method  that  requires  no  estimate  of 

\   .      or  \„.,r.   Since  \  .   and  \...„  may  be  arbitrarily  close  to  unity  (see 

mm     MAX         mm      MAX  * 

Lemma  k-3)   it  is  impossible  to  develop  a  general  a  priori  method  of  selecting 

an  initial  interval  (an  ,b,  )  c  (\   .  ,\,,.„)  which  avoids  the  inefficiency  of 

1'  1      mm'  MAX 

a  point  interval.  We  can,  however,  present  a  heuristic  method  of  selecting 

a  nondegenerate  initial  interval  that  yields  a  "good"  rate  of  convergence. 

The  interval  may  fail  to  satisfy  (a-,  ,b,  )  C  (\  .  ,X,,.__).   However,  it  will  in 

*   l  1     oan  MAX        ' 

this  case  yield  a  rate  of  convergence  defined  by  the  user  to  be  "good". 

Our  method  of  selection  is  based  on  the  following  facts.   First, 
Nnin  -1  —  \iAX*   An<3"  secon<^-y>   the  L  -norm  of  the  error  vector  is  reduced 

by  a  factor  o  <  — ; — -—- ,  after  N  iterations  when  (a,b)  ->  (\  .  .  \t.J)   is 

■"  T  (  a+b  ]  v  '  '        x  mm'  MAX^ 

Nl  b-a  ) 
used  as  a  constant  interval  in  iteration  (J*. 2).  We  will  show  that  if 
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(l)   a,  and  b  are  chosen  from  small  intervals  to  the  left  and 
right  of  unity  respectively. 


and 


ri-aJ 

where  the  user  decided  that  a  reduction  of  the  error  by  a 
factor  less  than  or  equal  to  €  is  "good", 

then  the  rate  of  convergence  will  be  either  the  "good"  rate  selected  by  the 

user  or  a  rate  which  is  optimal  or  nearly  optimal. 

There  are  three  cases: 

(i)  (a-,  ,b,  )  cz  (\   .   X.....)  , 
\  /   ^1*1'    v  min'  MAX  ' 

(ii)   (a,  ,bn  )  3  (\  .   *.,,.„)   and 

v  '      x   1'  1/    v  min'  MAX 

(iii)   eu,  <  \   .   and  b,  <  \,.AV  or  a,  >  \  .   and  bn  >  \„n„  . 

v    '   1    mm      1    MAX       1    mm      1    MAX 

If  (a,  ,bn  )    cz  (x   .    ,   \rir.„)   then  the  intervals  are  improved  by  the  methods  of 
v  1'  ly        v  min'   MAX 

the  previous  section,  and  approach  the  optimal  interval.   Thus  the  rate  of 

convergence  approaches  its  optimal  value.   It  should  be  noted  that  since 

(a-.  ,b..  )  c   (x   .        X...-U-)  we  will  have 
VT.'  ly        v  min'     MAX' 


< 


Tjal+bll  WNnin^MAX 


I  bn  -a,  /  I  X^.-^-K    .       I 

\    1     1 1  \    MAX     min  / 

Thus  the  optimal  rate  of  convergence  for  this  problem  is  slower  than  the  "good" 

rate  selected  by  the  programmer.   Also  note  that  since  X   .       Z   1  <  X.,,.,,, 
J  *     &  min  -   -  MAX' 

selecting  a,  and  b  from  small  intervals  to  the  left  and  right  of  unity  is 

a  good  heuristic  for  satisfying  (a,  ,b_,  )  a  (\  Xr,.v).      Comments  on  exact 

J   to  v  1'  1/     min'  MAX 

numbers  to  use  as  a,  and  b  are  given  below. 
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If  fa  i"b, )  =3  (\  .    ,   X„Av)  then  the  method  never  changes  the 
v  1'  1     min   MAX 

interval  and  the  rate  of  convergence  is  the  "good"  rate  selected  by  the 
programmer. 

Finally,  if  a^X^  and  b1>\m  or  a^X^  and  \<Xmx   then 
the  endpoint  in  (X  .  ,  X  hy)   is  improved  to  its  optimal  value.   The  endpoint 

ir0..ll    lYlriA 

not  in  (X  .    ,  X.„,J  remains  unchanged.  As  in  the  case  in  which  both  a,  and 
min   MAX  1 

b  are  in  (X   .    ,  X,  ,»„,.)  the  rate  of  convergence  is  improved.   This  improved 
1         min7  MAX 

rate  approaches  a  "nearly"  optimal  rate.   To  make  this  more  precise  consider 

the  case  a,  >  X   .      and  b  >\* Ay   The  left  endpoints,  a.,  will  be  decreased 

to  approximately  X   .  .   The  right  endpoint s,  however ,  will  never  be  improved. 

The  rate  of  convergence  will  approach  the  rate  obtained  using  the  constant 

interval  (X  .    ,  b_ ).   But  b,  was  chosen  in  a  small  interval  to  the  right  of 
mm   11 

unity,  say  b1  =  1  +  5.   In  addition  X^  >1  so  that  |  ^  -^  =  h±  -X^^,  <6. 
Thus,  intuitively,  although  b,  >X^        b,  is  not  a  gross  overestimate 
of  X         and  convergence  is  not  significantly  slower  than  the  optimal  rate. 

Remarks :  The  values  of  a,  and  b,  are  selected  to  yield  a  "good"  rate  of 

convergence  as  a  precaution  for  cases  in  which  (a,  ,b.  )  d  (X   .  ,  Xr,.v) .      The 

v  1'  1'  T  v  min-'  MAX' 

only  purpose  in  selecting  this  "good"  rate  is  to  prevent  the  selection  of 
too  large  an  interval  with  an  unacceptable  rate. 

It  is  impossible  to  create  a  rate  of  convergence  that  is  better 
than  the  optimal  rate  of  the  problem.   That  is,  selecting  (a, ,b_ )  such  that 


aN(a1,b1)  =     * 


VaJ 


is  very  small  cannot  improve  the  rate  of  convergence.   In  fact,  as 
o'N(a1,b1)  goes  to  zero  the  interval  becomes  a  point,  which  is  the  exact 
condition  we  are  attempting  to  avoid. 
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Our  experience  is  that  a.,  e  (l/k,   8/l0)  and  b   e  (1.25,  3)  yield 

"reasonable"  initial  intervals.   The  interval  (1/3,  5/2)  has  proved  to 

be  effective  for  problems  in  which  nothing  is  known  about  X   .      or  X 

mm     MAX 

In  Chapter  5  a  value  for  a,  in  some  special  cases  is  derived.   In  these 
cases  the  interval  (a,,  5/2)  has  been  effective. 

The  following  lemma  establishes  the  basis  for  the  heuristic  method 
of  selecting  the  initial  interval. 

Lemma  1+ .  3 :   Let  A  be  as  in  (2.6)  and  let  B  be  given  by  (2.7).   Then  1  is 
an  eigenvalue  of  (A+B)~  A.   Furthermore  for  every  e  >  0  there  is  some  A 

such  that  1-X  .  ((A+B)-^)  <  e  and  \MAV(  (A+Br^)  -1  <  e  . 

min  MAX 

Proof:   If  B  is  given  by  (2.7)  then  the  first  row  and  column  are  null. 

Hence  B  has  a  zero  eigenvalue.   Let  w  be  the  corresponding  eigenvector. 

Then  (A+B)w  =  Aw,  i.e.,  (A+B)~  Aw  ~  w. 

To  see  that  X   .      and  X....r   can  be  arbitrarily  close  to  unity 
min      MAX  J  J 

consider  the  following.   As  the  elements  B.   in  matrix  A  (2. If)  tend  to 

J  >*■ 

zero,  elements  b.   and  f    in  L  and  U  (2.6)  also  tend  to  zero.   In  the 

limiting  case  B.  ,  =  0,  b    =0  and  f .  ,  =  0.   Moreover  LU  =  A,  that  is, 

B  is  the  null  matrix.   In  this  case  (A+B)~  A  =  A~  A  =1.   Hence  every 

eigenvalue  of  (A+B)  A  is  unity.   Since  the  eigenvalues  vary  continuously 

with  the  elements  of  the  matrix,  for  a  non-null  B  with  arbitrarily  small 

elements  X   .      and  \.r.v   may  be  arbitrarily  close  to  unity, 
min      MAX 


Remark :   This  lemma  is  valid  for  the  factorizations  due  to  Kim  [11]  and 
Bupont,  Kendall  and  Rachford  [5]. 

We  shall  now  give  a  statement  of  the  algorithm. 
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The  Algorithm. 

1.  Select  an  initial  interval  (sUjlO  according  to  one  of  the  following 
rules. 

(a)  Define  (a..  ,b  )  according  to  the  method  presented  in  the 

first  part  of  this  section.  For  example  a,  €  (.25,  .75)  . 
and  b.  e  (1.5;  3.0)-  In  practice,  a1  =  l/j5,  b_  =  5/2  have 
been  effective. 

(t>)  If  possible  use  the  approximation  of  a,  developed  in 
Chapter  5,  Theorem  5-3*   Select  b..  as  in  (a). 

(c)  Let  a,  =  b  =  1. 

Comment :  Although  (c)  yields  a  convergent  and  practical- 
it  eration,  it  is  recommended  that  (a)  or  (b)  be  used  to 
define  the  initial  interval. 

2.  Select  an  initial  vector  x0  and  some  e>0  which  is  a  bound  on  the 
maximum  acceptable  error,  measured  by  the  norm  of  the  residual. 

3.  Select  an  integer  N  which  is  the  number  of  iterations  to  perform 

before  attempting  to  correct  the  current  interval. 

Comment:  As  N  is  increased  the  convergence  of  a.  and  b.  to  \  . 
11     min 

and  ^-MAy  is  improved  but  the  convergence  of  x  to  x  is  retarded  if 

the  current  values  of  a.  and  b.  are  crude  approximations  to  \   . 

11  ^  min 

and  a.^..      In  experiments  using  the  algorithm,   values  of 
N  =  5,  6,  7,  8,   9,   and  10  were  equally  effective. 
k.     Execute  N  steps  of  iteration  (k-2)  generating 

Ai       /\i  ,       Ai 

XN   '    8N-1      and      ^-1    ' 

where  i  corresponds  to   (a.,b.). 
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5.  Iflig 


All 


q-Ax  II  <  e  then  stop,  otherwise  compute 


11, 


in. 


1    <^*K.i$-i> 


u.   = 


6.    i.  If  ii.   <  a. ,      define   a  .    _    =  min  (min  a.  ,   /i.  )  and  b  =  b. 


i+1 


i+1 


j<i  y    x 

Add  one  to  i  and  go  to  k. 

If  jLt .  >  h.  ,  do  one  iteration  of  (k.l)   with  x  =  l/ju.  .   Define 

a.  ,  =  b.  and  "b.  ,  =  ju.  .   Add  one  to  i  and  go  to  k. 
l+l    i      l+l    i  ° 

Comment:   The  definition  of  the  next  interval  when  u..    >  b. 
i    i 

is  different  from  the  definition  when  ju.  <  a.  for  the 

l    i 

following  reason.   Components  of  the  error  vector 
corresponding  to  eigenvalues  X,    >  b.  will  have  been 


amplified  if 


3X> 


>  1,   where  the  i  indicates 


correspondence  to  the  interval  (a.,b.).   When  \i.    >  b. , 

the  next  interval  is  defined  to  emphasize  the  annihilation 

of  the  dominant  components  of  the  error.   The  components 

which  became  dominant  when  ju.  <  a.  ,  on  the  other  hand,  can 

ii 

become   dominant  only  if  they  are  not  diminished  as  rapidly 
as  other  components.      They  are  never  amplified  since 


PH(x) 


<  1  for  all  x  e  (0,a.).   Therefore  when  p..   <   a. 

'  l  ii 


it  is  unnecessary  to  concentrate  on  the  reduction  of  the 

dominant  components. 

If  u.    e  [a. .b. ]  ,  do  one  iteration  with  t  =  l/b.  .   Add  one  to 
i     i'  l   '  '  l 

i  and  do  N  iterations  according  to  ('4.22),  then  go  to  5- 
Comment :   This  requires  essentially  no  new  code  since  (>4.22) 
is  identical  to  (1+.2)  except  that  the  subscripts  on  the 
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right  side  of  the  equations  are  decreased.  Furthermore 
the  code  required  to  form  (A+B)_1(q-A^)  and  (i,  t(A+B)~  A)a£J 
before  proceeding  with  (4.22)  already  exists  in  the  code 
to  perform  (k.2). 

The  purpose  of  this  step  is  to  correct  the  condition 
as  described  preceeding  Theorem  if. 8. 


vw 


PN^MAX^ 


Convergence :  The  convergence  of  the  algorithm  is  established  by- 
considering  three  types  of  initial  intervals  (a,,b-). 

1.   If  (a.,  ,b,  )  ^  Ck  .    ,   \,,.,r)  then  the  algorithm  reduces  to  the  method 
v  1'  1'     mm   mAX 

of  Dupont,  Kendall  and  Rachford  E?]  and  is  therefore  convergent. 


2.   If  (a,,bj  i>  (X   .    ,   \.,AV)  and  P.A  .  )  =   P_T(\.AV) 


then  Theorem  If.  5 


guarantees  the  convergence  of  the  intervals  (a. ,b. )  to  an  interval 


1'  1 


(a.b)  such  that  a  =X  .   or  b  =  X.,.v.   Furthermore  if  (a,  ,K')c(\  .  A,,.,,) 
v  '  '  mm        MAX  1  1     mm'  MAX' 

then  the  (a.  ,b.  )  approach  (\  .  ,  X,,.,,)  and  the  rate  of  convergence 
^  V  \       *  mm'  MAX 

of  x  to  x  approaches  the  optimal  rate. 


3.   The  case  P,_(\  .  )  = 

x  N  mm' 

0  <  P  (\  .  )  <  1. 
Nv  min7 


PN^MAX^ 


is  trivially  disposed  of  since 


Work:  The  following  considerations  show  that  if  N"  is  the  number  defined 
in  step  3  of  the  algorithm  the  amount  of  additional  computation  per  iteration 
required  by  the  new  method  as  compared  to  the  method  of  Dupont,  Kendall  and 
Rachford,  (if. 2),   is  approximately  y^  times  the  number  of  computations 
in  one  iteration  of  their  method. 

There  are  two  differences  in  computation  of  the  methods.   The 

first  is  that  our  algorithm  computes  an  approximate  eigenvalue  on  every 

th 
N   iteration.   This  requires  2  vector  additions  and  the  accumulation  of  2 
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inner  products.   Therefore  the  number  of  extra  calculations  required  per 
iteration  to  compute  the  approximate  eigenvalue  is  (3k  additions  +  2K 
multiplications )/N,  where  K  is  the  dimension  of  A. 

The  second  difference  is  that  after  every  N   iteration  the  new 
algorithm  employs  a  new  interval.   This  creates  a  difference  in  the  work 
performed  because  the  first  iteration  with  each  interval  is  a  special  case. 
In  general 


4 


T  (y)        ,       T  ,(y) 


where  §  =  (A+B)~  £  ,  and  y  =  (a.+b.  )/(b. -a.  ).   When  n=0,  however, 
n  n-1  1  1  '   1  1  '  ' 


1  1 

is  computed.   In  every  case  x"  ,  is  defined  as  x"  +  S&  .   Thus  the  first 
*  n+1  n     n 

iteration  with  each  new  interval  requires  K  fewer  additions  and  K  fewer 
multiplications.   The  additional  calculations  due  to  these  differences  is 

(2K  additions  +  K  multiplications )/N 
The  number  of  computations  in  one  iteration  of  (J+.l)  or  (if.  2) 
without  the  calculation  of  the  approximate  eigenvalue  is 

(I3K  additions  +  UK  multiplications). 
Considering  only  the  multiplications  we  see  that 


the  additional  calculations 
per  iteration  of  the  new 
algorithm 


1 
UN 


calculations  per  iteration 
of  0.1)  or  (if. 2) 


if.  6  Experimental  Results 

Four  distinct  problems  will  be  studied  in  this  section.   Stone 
[l6]  and  Taranto  [17  ]  have  also  made  studies  using  the  type  of  problems 
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considered  here.   They  compared  the  amount  of  work  required  to  solve 
the  system  of  equations  Ax=q  using  first  Stone's  factorizations  and 
then  using  AD I.   Taranto  also  compared  Stone's  symmetric  factorization 
to  Stone's  nonsymmetric  factorization  and  studied  the  effect  of  changing 
the  parameters  a  and  x    (see  Chapter  2  for  definitions  of  a   and  t). 

The  following  tests  make  two  comparisons.   First  the  number  of 
iterations  and  the  amount  of  work  required  to  solve  systems  of  equations 
using  Stone's  symmetric  factorization  and  the  algorithm  developed  in  this 
chapter  is  compared  with  the  number  of  iterations  and  the  work  required 
when  the  problems  are  solved  using  Stone's  factorization  with  a  single 
best  parameter.   Next  the  comparisons  are  made  for  problems  solved  by 
Stone's  symmetric  factorization  with  the  new  algorithm  and  both  of  Stone's 
factorizations  used  with  Stone's  algorithm  to  vary  the  auxiliary  matrix, 
[15,  16]. 

Problems  a  and  c  described  below  are  identical  to  problems 
studied  by  Stone  and  Taranto.   Problem  b  is  similar  to  a  problem  they 
studied  but  differs  in  the  ratio  of  the  functions  a, (x)  and  a~(x). 
(see  Chapter  2,   where  a,  and  a.     are  defined.  )  This  change  was  made 
because  the  problem  studied  by  Stone  and  Taranto  converges  so  rapidly 
using  the  methods  under  consideration  that  no  significant  difference 
in  the  relative  rates  of  convergence  results.   The  last  problem 
studied  involves  random  numbers  in  the  matrix  A.   Since  the  random 
numbers  generated  in  the  problem  below  are  not  the  same  as  those 
generated  in  the  previous  studies,  no  comparison  of  results  is  attempted. 
In  fact  the  form  of  the  problem  is  changed  slightly  since  the  comparison 
is  impossible.   In  the  problem  studied  below  the  functions  a-,  and 
a2  of  (2.1)  are  random  throughout  the  region.   In  the  problem  studied 
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by  Stone  and  Taranto  the  functions  were  constant  on  subregions  but  were 
otherwise  random. 

Each  problem  is  solved  on  the  unit  square,  {(x, ,xp)  :  0<  x.  <  1  i=l,2), 
with  zero  boundary  condition,  i.e.,  Y(x)  of  (2.1)  is  identically  zero.   The 
solution  x  is  chosen  to  be 

x.   =  cos(jHh).cos(kHh)   j,k=l,  ...,  30, 

where  h  =1/(30+1).   Similar  tests  were  run  with  cosine  replaced  by  sine.   The 
results  were  almost  identical  and  therefore  are  not  presented.   The  vector 
q  of  (2.2)  is  computed  from 


The  initial  vector,  xn,  is  chosen  to  be  the  zero  vector. 


Problem  a.   This  problem  is  the  model  problem,  defined  by  letting  a, (x) 
and  ap(x)  be  identically  -1.   The 
by  the  analysis  of  Chapter  5  is  \. 


and  a0(x)  be  identically  -1.   The  initial  estimate  of  X   .  ,  suggested 
2  rain 


Problem  b.   For  this  problem,  called  the  generalized  model  problem, 

a,  (x)  =  -I/9  and  a~(x)  =  -1.   The  initial  estimate  of  X   .      as  suggested 
1       '       2  mm 

in  Chapter  5  is  5/8 • 

Problem  c.   This  problem  is  defined  by  a  heterogeneous  region  containing 
homogeneous  subregions  see  Figure  1+.2.   In  region  A,  the  ratio  of  a  (x)  to 
ap(x)  is  unity  with  a, (x)  =  -1.   In  region  B,  this  ratio  is  10  with  a, (x)  =  -1; 
in  C  it  is  0.1  with  a,  (x)  =  -l/lO  and  in  D  it  is  1  with  a..(x)  =  -l/lO.   The 
theory  of  Chapter  5  does  not  yield  an  estimate  of  X   .      in  this  case. 


6k 

Problem  d.   In  this  problem  random  values  of  a.  (x)  are  chosen  from  the 
set  {z  :  0<z<l)  . 

Each  of  the  four  problems  is  solved  several  times  with  N  of 
step  3  of  the  algorithm  equal  to  6  and  various  numbers  selected  as  initial 
endpoints  a  and  b  .   The  set  of  endpoints  was  chosen  so  that  the 
performance  of  the  algorithm  could  be  examined  for  both  good  and  poor 
values.   The  endpoints  were  chosen  to  create  the  following  situations. 

1.  a,  a  poor  approximation  to  X   .    ; 

2.  b,  a  poor  approximation  to  K^.yl 

3.  a,  and  b  poor  approximations  to  X    .      and  \»,AY  respectively. 
If.   a  and  b  chosen  according  to  the  suggestions  in  step  1 

parts  a  and  b  of  the  algorithm 
5.   Values  of  a  between  the  poor  and  the  suggested  values 
with  b.  =  2.5. 

The  problems  were  then  solved  using  a  single  best  value  of  t. 

Theoretically  this  value  of  t  is  2/(x   .     +KK.V);   however,  this  value  was 

'   mm   MAX  ' 

determined  experimentally.   Values  of  X   .   and  T^,.,,.  were  not  computed 

mm      MAX 

directly. 

The  results  of  these  tests  are  displayed  in  the  following 
graphs.   Only  linear  approximations  to  the  curves  are  shown  in  order  to 
simplify  comparisons. 


'  PN  (  XMAX  ) 


(b) 


/ 


PN(  ^MAX  ) 


^MAX5 


(C) 


Figure  l|.l. 

(a)  Two  cases:   In  Both  Cases  the  Component  Corresponding  to  \,«y  Is 
Dominant.   In  One  Case  the  Error  Is  Increasing,  in  the 
Other  It  Is  Decreasing. 

^   ^N^min^'^N^MAX^'  and  the  Err°r  Is  Decreasing* 
(c)  The  Component  Corresponding  to  X   .      Is  Dominant  and  the  Error 
Is  Decreasing. 
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Figure  4.2.  Heterogeneous  Test  Problem  Geometry 
(Region  Defined  on  Unit  Square) 
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Figure  4.3.   Model  Problem 
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Figure  4.4.  Generalized  Model  Problem 
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Figure  4.5.  Heterogeneous  With  Homogeneous  Subregions 
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Figure  1+.6.  Random  Numbers  in  the  Matrix  A 
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The  next  four  graphs  display  the  results  of  problems  a  and  c 
above  when  solved  by  the  following  three  methods. 


1. 


2. 


The  algorithm  of  Chapter  k   with  the  initial  interval 
(1/3,  5/2). 

Stone's  nonsymmetric  factorization  implemented  with 
a  sequence  of  a's  as  suggested  in  [l6],  namely 


m 


°kAX 


2h 


2h 


1+  a2/  1     1+  1/2 


(^. 31) 


and  a,  =  a  .   =  0:  the  intermediate  parameters 
1    mm    ' 

a.  ,   1  =  2,    . .  .  ,   m-1  are  computed  from 

i-1 
a     =   1-fl-a   )  m_l 


a.  and  ap  in  (U.31)  are  the  constant  values  of  a, (x) 
and  ap(x)  respectively  (see  Chapter  2  for  definition 
of  a  (x)  and  ap(x)).   In  problem  c  there  is  a  value 
of  a.    defined  for  each  subregion  using  the  local 
constant  values  of  a, (x)  and  ap(x).   For  the  tests 
presented  here  m=9-   The  a.    are  used  cyclicly  with 
each  a.    used  on  two  consequtive  iterations.   On  the 
odd  numbered  iterations  the  gridpoints  are  ordered  in 
the  normal  ordering;  i.e._,  the  gridpoint  (i,j)  (see 
Chapter  2)  precedes  the  gridpoint  (k,j)  if  i  <k  and 
precedes  the  gridpoint  (lj,m)  if  j  <m.   On  the  even 
numbered  iterations  the  reverse  ordering  is  used;  i.e._, 
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the  gridpoint  (i,j)  precedes  the  gridpoint  (is., 3) 
if  i  <k  and  precedes  the  gridpoint  (l,m)  if  j  >m. 
3.   Stone's  symmetric  factorization  implemented  with  a 
sequence  of  a' s  as  in  2. 

The  results  for  methods  2  and  3  are  taken  from  the  study  done 
by  Taranto  [17]. 

The  graphs  in  Figures  4. 7  and  4.8  have  the.  reduction  of  the  error 
is  plotted  against  the  number  of  iterations.   The  reduction  of  the  error 
is  calculated  as 

l-0«  /  HeJI 

where  e~  is  the  error  in  the  initial  vector  and  e  is  the  error  in  the 
0  n 

n   iterate.   The  graphs  in  Figures  4-9  and  4.10  have  the  reduction  of  the 
error  plotted  against  the  work  performed.   One  unit  of  work  is  taken  to 

he  the  amount  required  to  do  one  iteration  of  (4.2).   It  is  assumed  that 

cc 
all  the  auxiliary  matrices,  B  _,  are  calculated  only  once  and  that  the 

amount  of  work  needed  to  calculate  the  elements  of  the  factor  L  and  Tj  is 

equivalent  to  the  amount  required  to  do  one  iteration. 
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Figure  4.7.  Model  Problem 
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Figure  4.8.   Heterogeneous  with  Homogeneous  Subreglons 
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Figure  4.9.   Model  Problem 
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Figure  4.10.   Heterogeneous  with  Homogeneous  Subregions 


77 
Conclusions :   The  tests  comparing  the  algorithm  developed  in  this 
chapter  to  the  use  of  a  single  best  f  clearly  indicate  that  fewer 
iterations  are  required  when  the  new  algorithm  is  used.   Since  the  tests 
were  performed  with  N  of  step  3  in  the  algorithm  equal  to  6,   the  work 
estimate  of  the  last  section  implies  that  the  additional  work  required 
by  the  new  algorithm  is  approximately  l/66  of  the  work  performed  by  the 
iteration  which  doesn't  change  the  interval.   Thus  the  same  graphs  may  be 
used  to  conclude  that  the  new  algorithm  requires  less  work  than  the 
iteration  using  a  single  best  parameter. 

The  graphs  of  Figures  k>l   and  1+.8  indicate  that  in  some  cases 
the  required  number  of  iterations  of  the  algorithm  developed  is  greater 
than  the  number  required  by  Stone's  algorithm.   However,  Figures  k>9     and 
J+.IO  indicate  less  work  is  required  by  the  new  algorithm. 


5.  AN  ESTIMATE  OF  \    . 

mm 


Throughout  this  chapter  the  following  notation  will  be  used. 
B  will  denote  the  auxiliary  matrix  defined  using  (2.6)  and  (2.7)  with 
a   .  1.   The  quantities  ^  =  ^((A+B^A)  ,  b^fe  ,  c.^  ,  dj%k  , 

e.  ,  ,  f  ,  ,  C.  .  and  G.  .  will  be  as  in  (2.7).   B.  .  ,   D.  .  and 
j,k  '  j,k  '  j,k     j,k  j,k    j,k 

E    will  be  as  in  (2. If).   When  B.  .  ,   D    and  E.  ,  have  constant 
j,k  J^-K     j,  K      J,.K 

values,  except  B.   =  0  and  D.  ,  =  0  ,   they  will  be  denoted  by  B,  D 

and  E  respectively.   Whether  B  is  the  auxiliary  matrix  of  (2.6)  or  the 

scalar  constant  B.   of  (2. If)  will  be  obvious  from  the  use.   Finally, 
J  >k 

A 

the  elements  of  the  matrix  B  will  be  referred  to  as  b.  .to  avoid  any 
confusion  with  b.  ,  of  (2.6). 

A 

5.1  The  Limiting  Values  Of  b.  . 

iii 

In  this  section  it  will  be  shown  that  if  B.  ,  =  B  and  D.  ,  =  D, 

J,k         j,k 

A 

except  for  D,  ,  =  0  and  B    =  0  ,   then  the  elements  b.  .  ,   of  B,  approach 
J-^k        J?1  ijj 

a  limit  as  both  subscripts  increase.   In  particular  G .  ,  -*  G  =\IED   as  j  and  k 

J,k 

increase.   This  will  be  used  in  the  next  section  to  derive  an  estimate  of 

X    .    . 
mm 

The  following  lemma  will  establish  some  useful  properties  of 
A+B.  Although  we  are  only  concerned  with  the  case  a  =  1,  the  lemma  is 
valid  for  all  ae[0,l]. 


Lemma  5.1;   The  quantities  defined  in  (2.6),  (2.7)  exist  and  satisfy  the 
following  relations: 
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(a)  VSBj,k*° 

(b)  co,k*Vso 

(0)  dJ,k<0 

(d)  -1  s  e.^k  t=  0 

(e)  -1  *  f Jjfc  £  0 

(h>       cj,k  =  dj-l,k  ej-l,k 

j,k        j,k        j,k  "  j,k  ej,k-l  ~cj,k     j-l,k  £ 

Proof;    Parts  (a)  through  (h)  were  proved  by  Bracha  [1].   To  prove  (i) 
we  have, 


j,k    j,k    j,k  "  j,k  ej,k-l  "  cj,k  J-l,k 

j,k    j,kJ"  j,k  ej,k-l  "  Cj,k  j-l,k    j,k  "  j,k 

+  D    -  c.   -  b.  .  f .  .  n  -  c.  .  e 

J,k    j,k    j,k  j,k-l    j,k  j-l,k 

-  -V{1  +  ea,k.i  +  fj;k-i)  -  "j.kfr  +  ej.i,k  +  'j.i,k)  • 

Using  (a),  (b)  and  (f)  yields 

d.  ,  +  B.  ,  +  D.  ,  -  b.  .  e.  .  ,  -  c.  .  f .  n  .  §  0 

j,k    j,k    j,k    o,k  j,k-l    j,k  j-l,k 

which  completes  the  proof. 


Before  demonstrating  that  G  .  ,  =  c  .  .  f .  ,  .  tends  to  sT BD  as 

j,k    j,k  j-l,k 

j  and  k  increase,  it  is  shown  that  for  fixed  k  the  quantities  of  (2.7) 

approach  limits  as  j  increases. 
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Lemma  5-2:        Let  B.  .    =  B  and  D         =  D,    except  B         =  0  and  D  =  0   . 

r. J,.K  J,.K.  J, J-  -L,K 

For  fixed  k  the  elements  b.^  ,  c.^  ,  d.^  ,  e.^  ,  f  .^  ,  C.^  =  b.^  e.^ 

and  G    =  c    f  .  ,  ,  tend  to  limits  as  j  increases.   These  limits  will 
j,k    j,k   3-1, k 

be  called  bfc  ,  cfc  ,  d^.  ,  ek  ,   ffc  ,  Cfc  and  Gk  respectively. 


Proof:    The  proof  will  be  done  by  an  induction  on  k.   For  k  =  1  we  have 
the  following. 

b .  _  =   B .  .  =0  and  c .  _  =  D .  .  =  D  . 

3A   jji       oA   jA 


d 


n  =  _  2(B+D)-  De.  _  _  and  e.  _  =  D/d.  . 

hence 

d.A=-2(B+D)-  D2/d._1A  (5-l} 


)^I 


which  tends  to  a  limit  d  =  -(B+D)+Wb  +2BD 

Consequently  e  .  _  -*  en  =  D/d.  . 
3,1    1    '  1 

Finally  f    =  (B-Df.   -,)/d.   .   For  large  j  this  becomes 
3,1        J  ""1,1   3,1 

B-Df.  .  . 

f    _  J-1^ 

3,1        \ 

which  approaches  B/(d_+D)  as  j  increases.   Therefore  G.   -»  G  =  c  f 

Now  assume  the  truth  of  the  lemma  for  k  <  k.   For  k  =  K  we 
have  the  following. 


b  .  „  ->  B-G^  .  =  b^  . 
3,K      K-l    K 

c  .  Jr   -  D-G  .  =  c   . 
3,K      K-l    K 

d.  jr   -»  -2(B+D)+2Gt^  .  -  b  f  .  +  c  2/d   . 
3,K  K-l    K  K-l    K  '    K 

This  is  in  the  same  form  as  (5.1)  with  only  the  constants  changed.   Thus 

d,  „  tends  to 
3,£ 


81 


a  K"  -(MH-i-  T  \  <U  +  T"V(-2(BtD)+2\-i  "  \  fK-i'£-  **  2 


andej,K  ■*%  "  £ • 

Finally  fjj  -  %  -  B/(dK+  %)  and  8^  -  <fc  -  <fc  £  which 

completes  the  proof. 

The  next  lemma  will  show  that  the  quantities  b  and  cl  of  the 

previous  lemma  are  increasing  in  magnitude  as  k  increases  if  0  >  B  >  D. 
Lemma  5-k   then  shows  that  the  b,  are  increasing  in  magnitude  for  any 

B,  D  <  0.   This  result  is  used  in  Lemma  5-5  to  prove  that  as  j  and  k 
increase  the  quantities  t> .  ,  and  d.   are  approaching  limits  h  and  d 

respectively. 


Lemma 


5.3:   Suppose  0  >  B  >  D  .  Let  b  and  cl  be  as  in  Lemma  5-2.   Then 


the  b,  are  decreasing  and  the  d  are  increasing  as  k  increases. 


Proof;   The  proof  will  again  be  by  an  induction. 
b.  =  0. 


c   =  D. 


d   =  -(B+D)+\yB^+2BD  . 

f     =  B/(d     +  D)    =  (B  +VB2+2BD)/2D  . 

and 

G,  =  c  f   =  (B  +\/  B2+2BD)/2. 

To  show  that  b  <  b  and  d  >  d  we  have  the  following  relations. 


82 


b2    =  B  -    G„ 


C2 


=   (B^\jB2+2BD)/2  <  0    =  1d1     ? 
o  D  -  G     =  (2D-B'^\/B2+2bD)/2     , 


and 


d2    =  -(B+D)+G1-  -|>  b2   fx  +  ~  {[-2(B+D)+2G1-b2f1f   -l+c^}  2 
Substituting  the   above   expressions  for  f      ,   G      ,   b     and  C     yields 

d2    =—--]>  -|^\/b2+2BD  +  -|-   (-  -2-  B2  +  6BD-3B~\/B2+2BD)"2"     . 

Since  B   >  D, 

-3ByB2+2BD  >  -3b\/B2+2BD  +  ^p-  +   B2  +   2BD-6BD.  (5.2) 

Collecting  terms  taking  the   square  root  and  multiplying  by  l/2  we  obtain 
-i(-  -|-  B2+6bD-3BVB2+2BD)1/2>     -  -|-  B+|"WB2+2BD   . 

Using  this  inequality  and  the  last  expression  for  d    yields 
d2  >  -B-I>Vb2+2BD   =  d     . 

The  following  relations,   which  come  from  substituting  the  limiting 

values  of  b  and  d.        into    (2.7)   will  be  required  in  the  general  step  of 

3  )*■  J  )K- 


the   induction. 


b.    ,  -B+D     -1 
k-1 


\=B(1  +  \~T]  *=2>3>'"  (5-3) 

and  p 

b  (b   -B+D)^ 

Note  that  b  <  B  <  0  implies 
k 

b   -B+D 
and  thus 

-  k  4^ —  <  0.  (5.5) 


dk. 


1 
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Before  starting  the  inductive  step,  it  is  necessary  to  show 

that  b,  <  b0.   This  inequality  holds  if 
3    2 

or  equivalently  if 


b  -B-D 

<  "411-  •  (5.6) 


d2         dl 


Replacing  b^  and  d  by  (B^\/b  +2BD)/2  and  -(B+D)  +WB  +2BD  (5.6)  becomes 


-3BD-2D2+3D  B2+2BD 
2d^ <  -B  +  D  . 

Substituting  for  d  and  collecting  terms  we  get 


(B+2D)  (-  -|-  ^yB2+2BD)-  -i  BD  <  (D-B)  (-  -2-  B2+6BJ)-jS\J  B2+2BI))1'2 .     (5-7) 


This  will  be  true  if 

B2 


(B2+l+B]>lfD2)    (  -2-     -ByB2+2BD     +   B2+2BD)    -  7BD(B+2D)*Vb2+2BD 


is  greater  than 


(D2-2BB+B2)    (-  -£■  B2  +6BD-3ByB2+2BD)    . 


4- 

(That  D-B  is  negative  was  used  in  obtaining  this  condition. )  Multiplying 
out  the  factors  of  the  above  inequality  leads  to 

2B^-l/2  B5D  +25-2-  B2D2  +2BD5  -ByB2+2bd  (-2B2+15D2+17BD)  >  0 

which  is  obviously  true.   Thus  (5-7)  is  valid  and  therefore  (5. 6)  is  valid 

and  b^  <  b„  . 
3    ^ 

We  now  start  the  inductive  step.  Assume  b  <  b    and 


8k 


\-l  >  \-2  ■   ^  \   =  Vl  "  \  >  °  '   aM  let  ek  =  ^  "^-1  •  From  (5.U) 


b  2      (b.  -B+D)2 
k      v  k 


Replacing  bfc  with  bfe_1  -  ^  and  d^  with  €k  +  d^  introduces  -  2(l>f"bk_1)' 
Then  using  (5.0  and  the  definition  of  efc  yields 


Since  d.  2  <  d,  ,  by  assumption, 

^_;L+2  Wl'^  * B2+I)g-2Bbk_1+2I>bk_1-2BJ   (yyM)8 

e*>28* \7i \.x  ^ 

2&  e  (b  -B+D)2 

■^  K-A-*^' +      a^    •  (5-9) 

Now  from  part  (i)  of  Lemma  5*1  we  have  d,  ,+B+D-b,  ,  e  p-c,  ,f .  1  >  0. 

But 

-Viek_i  =  "ck-2fk-2  =  \-rB   > 


and 


-c.    _f.    _  =  b  -B 
k-1  k-1       k 


Hence 

-1  ■""""  "k-l~~'"k 


d  i: 


and 

t^VVi  >0  •  (5-10) 

Now  (5.9)  and  (5. 10)  together  with  5,  >  0  and  d.  ,  >  0  yield 

e  (b  -B+D)2 

k  >    K 


V4 


1  k 
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Furthermore  (5-10)  implies  |  (bk-B+D)/dk_1|  <1  and  (5.5)  implies 

I  Ob. -BKD)/d.  I  <  1.   Hence  e  >  0  and  therefore  d>  d 

It  remains  necessary  to  show  that  b,  _  <  b  .   We  will  show  that 

e  <  e    ,   where  e  is  defined  as  in  Lemma  5.2  (This  in  turn  requires 

fk-l  <fk-2'  (5*12)  bel°W  )' 
From  (5- 3) 


Vi 

=  B 

1 

bn  -B+D 

-1+  X 

But 


b.  -B+D 
k 


\        \ 


=  e, 


Thus  b_T  m      B  (  -rr- ) 

k+1        1+e,  ' 
k 


Therefore  if  e.,  <  e,  ,  then 
k    k-1 


(5-11) 


bk+l  <  B/^1+ek-l)  =  \  ' 

Since  b,  <  b .  .  ,    (5. 11)  implies  e ,  _  <  e,  „.   Next  from  the 


definition  of  d    along  with 
3  )*■ 


-2B+C  ,f,  ,+b,  f,  ,  =  -2c,  =  -2d  e 


'k-1  k-1  k  k-1 


V 


we  arrive  at 
2 


V*A   =  -2B-2dkek-Vk-i  • 


Thus 


(l+ek)   =  (-SB-V^J/^  . 
But  <1  >   d.  .  and  b   <  b   .   If  in  addition 


f    <  f 
k-1    k-2 


(5-12) 
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Then  ^fw  <  -^-\_±\_2   and 

_   -2B-b,  ,f,  p  p 

(1+e/  <     k-X  X-2   -  (1+e^)2  .  (5.13) 

But  l+en  and  l+en  ,  are  positive  so  that  (5.13)  implies 
k       k-1 

1+e.  <  1+e  _ 
k      k-1 


or 


e,  <  e.  n  .  (5.H) 

k    k-1 

To  see  that  f.    .   <  f ,    n  we  write 
k-1  k-2 

k     k-1  k-1  k-1     k-2  k-2 

Since   c,  -e,    ,    =  b,  -b,    _    ,   the  last  equation  implies 
k     k-1         k     k-1 

•WWA-i  -  ck-i  +(°k.i+8k.i)fk-e  • 

"w'Vl^    "  5k+8k-lfk-2   • 

Vi(fk-rfk-2»+(-M)(fk.rfw»  -  VV-Ajs  • 

fk-2<-\-2+B-D'   >  fk.l(-Vl+M)    • 


But  c.  =bn  -B+D   so 
k     k 

-c     , 

f         >  f  £Si_  >  f 

k-2         k-1     -c,    _  k-1 

k-2 

Hence  (5-12)  holds  and  implies  (5-llj.)  which  implies  b    <  b  .   This 
completes  the  proof. 


If  we  remove  the  condition  B  >  D  we  can  still  prove  b    <  b 
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for  all  k  although  it  is  no  longer  true  that  d.  ,  >  d,  .   The  restriction 

B  >  D  was  used  only  in  (5-2)  to  show  d  >  d  and  in  (5.7)  to  show  b  <  b  . 

If  D  ^  B  then  (5*7)  is  trivially  true  and  thus  b,  <  b_  .   Furthermore  if 

for  some  k,  d,  ^  cL,  and  b  <  b    ,  then  from  (5.3) 


b  -B+D   \  _1   /    b,   -B+D  \   _1 

That  is,  b    <  bk  whether  d^^dkor<L>cL.   Thus  the  previous  lemma 
handled  the  more  difficult  case  of  the  following  lemma. 


Lemma  5>k:         The  b  as  defined  in  Lemma  5»2  are  decreasing. 

Lemma  5.U  can  now  be  used  to  prove  all  the  elements  of  the  matrix  B 
approach  limits. 


Lemma  5 . g :   Let  B    =  B  and  D.   =  D.   The  quantities  b    ,  c.   , 
3  >&  3  >*■  3  >&■         3  >x- 

d.  .  ,  e.  .  and  f .  .  as  defined  in  (2.1)   approach  limiting  values  as 
J,k    j,k      j,k  **  *° 

j  and  k  increase.   We  will  call  these  values  b,c,d,e,  and  f  respectively. 


Proof:   Since  all  elements  can  be  defined  in  terms  of  b .  ,  and  d.  .  , 
J,k      j,k 

we  will  prove  only  that  b .  ,  -»  b  and  d .  ,  -*•  d. 

J  ,&         3  >*■ 

By  Lemma  5.1+,  the  b  are  increasing  in  magnitude  as  k  increases. 
It  is  sufficient  to  show  that  the  b  are  bounded  below  for  then  they  must 
approach  a  limit,  b.   Moreover  the  d^  must  approach  a  limit  since 
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b.     .-B+D        , 

bk  =  B(1+  ^_)- 


Let  d   .      =  rain  {&.  }    and  let  b,    =  0   ,   b      =  b     and 
mm  kj  1  d         d 


b   =  b  ( Ssia ) 

k  d   .    +b\     ,-B+D 

mm     k-1 


Now  suppose  b     ^  b     then 

d   . 


k+1  d   .    +b?  -B+D 

mm     k 


d   . 
mm 


d   .    +b,  -B+D 
mm     k 


since  the  denominator  is  positive  and  increased.   The  fraction  is  greater 
than  one  so  that  adding  the  same  positive  number  to  the  numerator  and 
denominator  decreases  it.   Hence 


\    n  =  B       \  =  b,  n  . 

k+1      vv»      k+1 

Thus  b\  ^  bn  for  all  k.   Furthermore  the  bn  decrease  to 

k         k  k 


-d   .    +  B-D+V  (d   .    -B+D)2  +l^Bd   . 

limh     =    212 V      minn 21HL     . 

k-*  oo 

The  b  ,   therefore,  are  bounded  below  and  the  proof  is  complete. 


The  next  theorem  establishes  values  for  the  limits  of  b,  and 


Theorem  5.1:    If  b  and  d  are  as  in  Lemma  5.2  then 


\ 
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lim  b      =  B- VbD     and     lim  cL     =  -B-D+2  V"bD  . 
k  -*  oo  k  «*  oo 

Proof:        From  Lemma  5-5  ^  have  that  b,    and  cL    approach  limits     b  and  d 

respectively.      In  addition  these  limits  satisfy  (5.3)   and  (5.k)   since 
the  b     and  cL    satisfy  them.      Hence 

Bd J^     b-B+D  \  -1  .  , 

b   =     d+b-B+D  =B     1+— (5.15) 

and 

-1        o  2 

d  =  -2(D+b)    -  i  (b  +  (b-B+D)    )    . 

Therefore 

d2+2d(D+b)+b2+ (b-B+D)2   =  0   .  (5.16) 

From  (5.15) 

b2   m  B2(l+     2(b-B>D)        +      (b-B+D)2        fl      _ 

d 

2  2 
Using  (5. 16)  to  replace  (b-B+D)  /d  yields 

b2  =  B2  (i+  2(b-B+D)  _  x  _   2(D+b)   _  _bf_  )   -1 


d  d        ,2 

d 


2  2 

B  d 


-2Bd-b2 


Consequently 


k   P         ?  P 

-b  -b  (2Bd)-B  d  =  0. 


and 


d 


2  2   2  2 
,2    2Bd  -  VkB  d  -UB  d 
b  = — ■  -Bd  . 

Since  b  <  0  ,   the  last  equality  implies 

b  =  -V^id   .  (5- IT) 
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Substituting  (5-17)  into  (5.16)  yields 

d2+2d(D-V^Bd)-2Bd  -2"V-Bd  (-B+D)+(-B+D)2  =  0. 
Collecting  terms  we  get 

(d-B+D)2-2V^Bd~(a-B+D)  =  0. 
Hence 

-2b    =  2"V  -Bd   =  d-B+D  .  (5-18) 

Squaring  both   sides  of  the   second  equality  and  again     collecting  terms 


0   =  d2+B2+D  +2Bd+2dD-2BD. 


Therefore 


d 


+\        2        2  2        2 

-2B-2D  -  V  1|3  +4D  +8BD+8BD-^B  -W> 
d   =     2 

=  _B-D   -2  VBD 

Note  that  if  d  =  -B-D-2  VED  then  (5. 18)  implies  b  =  B+V  BD  .   But 

b  ^  b_  <  B. 
k 

Therefore  we  must  have 

d  =  -B-D+2  VBD 
and 

b  =  B-VBD  . 
This  completes  the  proof. 


Corollary  5. 1.1;    If  B.  .  =  B  and  D.  ,  =  D  then  G.  ,  of  (2.7)  tends  to 
0,k         J,k  j,k 


=V^ 


G  =  yBD  as  j  and  k  increase. 


Proof:   This  is  obvious  from  the  theorem. 
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5.2  An  Estimate  Of  X. 


rain 


In  this  section  an  estimate  of  X   .      to  be  used,  in  the  algorithm 

mm  to 

of  Chapter  k  is   developed.     We  "begin  "by  examining  the  case 

B.  ,    m  I  and  D.  .    a  1  except  B.        ■  D,   .     =0. 
J,-K  o,k  o,l         l,K: 


Then  some  more  general  cases  are  considered. 


-1, 


It  will  he  shown  that  for  B  =  D  =1,  X   .   ((A+B)~  A)  approaches 


mm 


l/2  from  ahove  as  the  meshsize,  h,  of  the  grid  used  in  deriving  (2.2) 

from  (2.1),  tends  to  zero.   We  begin  with  the  following  lemma  which  will 

be  used  to  show  X   .      ( (A+B)  A)s  l/2  for  all  h  >  0. 

mm  ' 


Lemma  ^.6:       Let  A  be  as  in  (2.k)   with  B.  ,  =   D    =1  except  B.   =  D   =  0 
3  >&         J^-K  J,l    1,£ 

j,k  =  \,2,'"  ,   n  .   Let  B  =  ($.  .)  be  the  auxiliary  matrix  derived  from 

^-  >  3 

(2.7)  using  a  =  1.      Then  for  all  x  £  0,     <Ax,x  >  >    <Bx,x  >   . 


Proof:  Bracha  has  shown  in   [1]  that  if  A  =  (a.     .)   is  an  N  by  N  matrix 

in  the  form  of   (2.4)   and  B  =  (b\     .)   is  given  by   (2.7)   then 


1M-J. 

<Ax,x>      =     J       -a_ 


i+1 


|2  ¥  I 

x.,n-x.  +       y  -a.    .        x.      -x. 

1+1     1  .*•  1,1+n     i+n     1 

1  1=1  ' 


N 

>  (a.    .+a.    .  ,-,+a.    .,    +a.    .    .+a.        , 

£*.  1,1     1,1+1     i,i+n     i-l,i     i-n,l" 

1=1  ' 


x. 


and 


<Bx,x> 


V-       A  .  ,2  V       A 

=     -      Z.       b.    .  ,-.      x.  ,    -x.  +      )        b.    .  , 

f*n       1,1+1    j    i+n     i|  f*p       1,1+n 


x. 
1 


N-l 


N-l 


V  b.    . ,_    x. ,  ,-x.  +      Y         b. 

A*    ,       1,1+1     1+1     i  A»  ,,     1,1+1 

i=n+l  '  '  i=n+l       ' 


x. 

1 


(5.19) 
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N-n-1 


k    *"> 


i+n 


x. , , -x. , 
-l+l     i+n 


2  ^Ji-1 


+        £         bi+l,i+n(lXi+lf 
i=l 


+  x. 


i+n 


■), 


where  N  =  n 

The  second  and  fourth  sums  may  be  rewritten  as 


N-n-1 


I     e, 


i=l 


i+l,i+n+l 


Xi+1 


(5-20) 


and 


N-^"1  a  2 

>     h.    .    ,  x. 
*-•     i+n,i+n+l 
i=l 


i+n 


(5-21) 


Now  from  the  definition  of  B 


i+l,i+n+l 


A 

l+l, i+n 


and 

A  A  A 

"b.        .        .      =  -  ah.        .    .    =  -  ah.    _    . 
i+n,i+n+l  i+n,i+l  i+l,i+n 

Hence  for  a  =  1,    (5.20)   and   (5.21)  become 


N-n-1 


i=l 


i+1, i+n 


x 


i+1 


2  N"^-1  a 

and  -      )         D.,..    . 

A-         l+l,  i+n 

i=l 


x. 


i+n 


respectively.   Therefore, 


N-l 


<(A-B)x,x>=   Y       -a.    .    _ 


X.  ,  -  -X. 

1+1     1 


2      r1  ft 

i=n+l        ' 


x       _x 
i+1     i 


N-n 


i=l       ^1+n 


x. ,    -x. 
i+n     i 


2        N~n  / 
+ 


ft     L'1+n 


x.      -x. 
i+n     i 


(5.23! 
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N 
+  5^     (a.    ,+a.    ,.,+a.    ,,    +a.        .   +a  )!x 


N-n-1 


i=l 


i+l,i+n 


l+l     i+n 


N-l 


f^         i,i+l     i,i+l' 

N-n  A 

+    T.      (-a.    .      +b.    .      ) 
£■*,    v     1,1+ n     1,1+rr 


N-n 

t 

ltd 
N 

I 

i=l 
N-n-1 


x       -x. 

l+l     l 


X.  ,     -X. 

i+n     i 


+  >        (a.    .+a.    ,  ,,+a.    .,    +a.    -.    .  +  a.        .) 

«-*         1,1     i,i+l     i,i+n     i-l,i     i-n,i 


x. 


A 


T       b-  -i  • 


i+1 


Xi+l"Xi+l 


(5.25) 


From  Lemma  5-5  and  Theorem  5-1  it  is  seen  that     b 


1;J 


<  1  and  from 


Lemma  5-1  it  is  seen  that  b. .,  .   >  0.      Hence  all  the  coefficients 

H-l, i+n 

in  (5-23)  are  positive  since  a.  .  =  -1  for  i  ^  j.   Therefore  <(A-B)x,x>  > 
proving  the  lemma. 

Theorem  $.2;    If  h,  A  and  B  are  as  in  Lemma  5-6,  X   .  ((A+B)"  A)  >  l/2 


mm 


and  X   .    ((A+B)_1A)   -*  l/2   as  h  -»  0. 

mm"  ' 


Proof:     From  Lemma  5-6  we  have 
<(A+B)x,x>    <2  <Ax,x>  . 

^„,     <Ax,x> .       <Ax,x> 

ThUS     <(A+B)x,x>         >      2~4x^ 


=  1/2   .      Using    (3.6)    yields 


A    .    ((A+B)-1A)  >   1/2 
mm  ' 


X 


9h 
To  establish  the  remainder  of  the  theorem  consider  the  vector 

i  \T    ,  ,,  ,      _/l  if  j+k=n  j  >   J  and  k  >  K   /E  OJ  % 

=  (x    •••x   )   such  that  x.  ,  =(~   ,,    .  (.5. 2k. ; 

k  1,1'    n,ny  j,k  10  otherwise 


where  J  and  K  are  chosen  so  large  that  the  elements  of  B  have  nearly 
constant  value.   That  is,  if  row  i  of  B  corresponds  to  the  point  (j,k) 
with  j  >  J  and  k  >  K  then  b±  ±_^  =   -1  ;  \>±   i_n+1s  1  I   1°±   ^j-  -1  5 

b   =  2  ;  b.  .  .=  -1  ;  b.  .   .=  1  and  b.  .   =  -1. 
ii    '   1,1+1     '   1,1+ n-1        1,1+ n 

<Ax,x>  =  ( n-K-1  )><.  and  <Bx,x>  =    (n-K-3)^+6. 


<Ax,x>  _  ( n-K-1 K 1 ft,  pe-\ 

<(A+B)x,x>       "       2(n-K-l)4-2        ~     I  2  '  ^'^ 

2"      (n-K-l) 

as  h  goes  to  0,  (i.e.,  as  n  goes  to  infinity)  the  quantity  in  (5.25)  tends 

to  l/2  from  the  above.   Using  (3-6)  again  completes  the  proof. 


The  argument  used  to  prove  the  second  half  of  the  above  theorem 
can  be  used  to  prove  a  similar  result  for  a  more  general  matrix  A.   Suppose 
A  is  derived  from  a  differential  operator  in  (2.1)  with  constant  a, (x) 

and  a0(x).   Then  A  will  have  constant  diagonals  a.  .   =  B.  a.  .  _=  D 
2  D       1,1-n       i,i-l 

and  a.  .  =  -2(B+D)  except  for  some  zero  elements  on  rows  corresponding 

to  boundary  points.   Furthermore  the  elements  of  the  auxiliary  matrix  B 

tend  to  the  following  values:  b.  .+  ->  -V  BD  :  b.  .+  ,  „  N  -A/BD  • 

1,1-n  '      i,i-(n-l) 

b       +n    -»  -V  BD  and  b.    .    ->  2VBD  .      Thus  if  X  is   defined  as  in   (5-21+)   then 
1  ^i-J-  1 ,1 

<Ax,x>  (n-K-1)    (-2)    (B+D) 

<(A+B)x,x>  (n-K-1)    (~2)    (B+D)+ (n-K-3)My^  )^^~ 


-2 (n-K-1)    (B+D) 

2 (n-K-1 J    (-B-D+yBD)    -2yBD" 
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-(B+P) 1 

-(B+D)  +2VbF      p+  2VBT 

(3.6)  and  the  above  yield  the  following  theorem. 


Theorem  5.3:   If  A  is  in  the  form  of  (2.6)  with  B.  .  =  B  and  D.  ,  =  D 

except  B .  .  =  0  and  D    =  0  and  if  B  is  the  auxiliary  matrix  defined  in 

(2.10)  using  a  =  1,  then 

lim^aA+B)-1^  t-     \m  ,  (5.26) 

where  h  is  the  meshsize  of  the  grid  used  in  the  discretization  of  the 
differential  problem. 


There  is  no  lemma  coresponding  to  Lemma  5*6  which  implies  the 

inequality  in  (5*26)  can  be  reversed.   That  is,  the  quantity  in  (5.26) 

is  not  lim  X   .  ((A+B)~  A).  However  this  value  has  proved  to  be  effective 
h->  0  mn 

as  the  initial  left  endpoint,  a,  ,   in  the  algorithm  of  Chapter  k-      This 

value  is  suggested. 

The  result  of  the  previous  theorem  can  be  generalized  to  a  martix 

A  which  is  block  tridiagonal  with  constant  diagonals  within  each  block.  The 

result  is  of  little  value  however  since  again  there  is  no  lemma  which  like 

Lemma  5-6  establises  a  lower  bound  for  the  eigenvalues.   Furthermore,  the 

value,  which  arises  from  this  modification  of  the  argument  leading  to 

Theorem  5«3>  is  so  crude  an  estimate  of  X   .  that  it  is  not  worth  calculating 
^'  mm 

for  use  in  the  algorithm. 

As  mentioned  in  Chapter  k   the  initial  left  endpoint,  a,  ,  to  use 
in  the  algorithm  should  be  calculated  according  to  Theorem  5-3  when  possible. 
If  the  Theorem  is  not  applicable  then  a  e[l/l,l/2]  is  sufficient. 
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